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

#
# NEON 64-bit and 128-bit SIMD ISA Definitions
#
# legend: v = ok, ! = missing, u = added, untested
    # Required:
    #  v active
    #  v bin_shl1  bin_shl2  bin_shr1  bin_shr2  bin_shrev
    #  v bits, ctype, t, v
    #  v countrec
    #  v includes
    #  v info
    #  v instr
    #  v isFixedPoint, isFloat
    #  v loadCont (uses sv loads, ARM does not have masked loads)
    #  v mul_cx, mul_cx_conj
    #  u reverse
    #  v splopts
    #  v storeCont (uses sv stores, ARM does not have masked stores)
    #  v svload
    #  v svstore
    #  v RCVIxJ2
    #  v dupload, duploadn
    #  ! hadd
    #  v swap_cx
    #  v vzero

    # Fixed point:
    #  fracbits
    #  saturatedArithmetic

    # Viterbi:
    #  ! interleavedmask, hmin, average (?), isSigned

Class(SIMD_NEON, SIMD_ISA, rec(
    file := "neon",
    commonIncludes := self >> [],
    active       := true,
    isFixedPoint := false,
    isFloat      := true,
    ctype        := "float32_t",
    stype        := "__attribute__ ((aligned(16))) float32_t",
    bits         := 32,
    useDeref     := false,
    splopts      := rec(precision := "single"),
    arrayDataModifier := "__attribute__ ((aligned(16)))",
    arrayBufModifier  := "static __attribute__ ((aligned(16)))",
    declareConstants := true,
    threeOps         := false,
    fma		     := true,
    vzero := self >> self.t.zero(),
    unparser := NEONUnparser,
    compileStrategy := self >> BaseIndicesCS
      :: When(self.fma, [DoFMA], [])
      :: [ MarkDefUse, #
           (c, opts) -> CopyPropagate(c, CopyFields(opts, rec(autoinline := true))),
           MarkDefUse, #
           (c, opts) -> CopyPropagate(c, CopyFields(opts, rec(autoinline := true))) ]
      :: When(self.threeOps, [DoThreeOp], [])
      :: [ Compile.declareVars,
           (c, opts) -> opts.vector.isa.fixProblems(c, opts),
           HashConsts ],

    loadCont := (self, n, y, yofs, x, xofs, xofs_align, opts) >> let(
	nn := _unwrap(n),
	yy := vtref(self.t, y, yofs),
	assign(yy, deref(nth(x, xofs).toPtr(self.t)))),

    storeCont := (self, n, y, yofs, yofs_align, x, xofs, opts) >> let(
	nn := _unwrap(n),
	xx := vtref(self.t, x, xofs),
	assign(deref(nth(y, yofs).toPtr(self.t)), xx)),

    storeContAcc := (self, n, y, yofs, yofs_align, x, xofs, opts) >> let(
	a  := _unwrap(yofs_align),
	nn := _unwrap(n),
	xx := vtref(self.t, x, xofs),
	t  := TempVec(TArray(xx.t.t, xx.t.size)),
	decl([t], chain(
	    self.loadCont(nn, t, 0, y, yofs, yofs_align, opts),
	    assign(vtref(self.t, t, 0), vtref(self.t, t, 0) + xx),
	    self.storeCont(nn, y, yofs, yofs_align, t, 0, opts)))),
));

Declare(NEON_HALF);

Class(NEON_HALF, SIMD_NEON, rec(
    info := "NEON 2 x floating point",
    v     := 2,
    t     := TVect(T_Real(32), 2),
    freshU  := self >> var.fresh_t("u", TVect(self.t, 2)),
    freshT  := self >> var.fresh_t("t", self.t),
    instr   := [vunpacklo_half,  vunpackhi_half, vrev_half ],

    includes := self >> self.commonIncludes() ::
        ["<arm_neon.h>", "<include/omega32_neon.h>", "<include/mm_malloc.h>"],

    dupload := (y, x) -> assign(y, vdup(x, 2)),
    duploadn := (y, x, n) -> assign(y, vdup_lane_half(x, n)),

    reverse := (y,x) -> assign(vref(y,0,2), vrev_half(vref(x,0,2))),
    RCVIxJ2 := (y,x,opts) -> assign(y, vrev_half(x)),

    # this is shift "right" in Spiral notation (and ARM too)
    # [a b] -> [0 a]
    bin_shl1 := (self, y,x,opts) >> assign(y, vext_half(self.t.zero(), x, 1)),
    # this is shift "left" in Spiral notation (and ARM too)
    # [a b] -> [b 0]
    bin_shr1 := (self, y,x,opts) >> assign(y, vext_half(x, self.t.zero(), 1)),

    # this is shift "right" in Spiral notation (and ARM too)
    # [a b], [c,d] -> [b c]
    bin_shl2 := (self, y,x,opts) >> assign(y, vext_half(x[1], x[2], 1)),
    # this is shift "left" in Spiral notation (and ARM too)
    # [a b], [c,d] -> [b c]
    bin_shr2 := (self, y,x,opts) >> assign(y, vext_half(x[1], x[2], 1)),

    # [a b], [c,d] -> [b c]
    bin_shrev := (self, y,x,opts) >> assign(y, vext_half(x[1], x[2], 1)),

    countrec := rec(
        ops := [
           [ add, sub ],
	   [ mul, vmulcx_half ],
           [ fma, fms, nfma ],
           [ vunpacklo_half, vunpackhi_half, vrev_half, vswapcx_half, vtrnq_half, vextract_half,
	     vdup_lane_half ],
           [ vload1_half ],
           [ deref, nth, ],
           Value      # Value without [] is a keyword in countOps !!
        ],
        printstrings := ["[adds]","[mults]","[fmas]","[vperms]","[svldst]","[vldst]","[vval]"],
        arithcost := (self, opcount) >> opcount[1]+opcount[2]
    ),

   loadCont := (self, n, y, yofs, x, xofs, xofs_align, opts) >> let(
	nn := _unwrap(n),
	yy := vtref(self.t, y, yofs),
	p0 := nth(x, xofs).toPtr(self.t.t),
	Cond(
	    nn=1, assign(yy, vload1_half(p0, self.vzero(), 1)),
	    nn=2, assign(yy, deref(nth(x, xofs).toPtr(self.t)))
	)
    ),

    storeCont := (self, n, y, yofs, yofs_align, x, xofs, opts) >> let(
	nn := _unwrap(n),
	xx := vtref(self.t, x, xofs),
	p0 := nth(y, yofs).toPtr(self.t.t),
	Cond(
	    nn=1, vstore1_half(p0, xx, 1),
	    nn=2, assign(deref(nth(y, yofs).toPtr(self.t)), xx)
	)
    ),

    svload := [
        # load using subvectors of length 1
	[   (y,x,opts) -> assign(y, vload1_half(x[1].toPtr(NEON_HALF.t.t),
		                                NEON_HALF.vzero(), 1)),

            (y,x,opts) -> let(u := var.fresh_t("T", NEON_HALF.t),
                decl([u], chain(
                        assign(u, vload1_half(x[1].toPtr(NEON_HALF.t.t),
				              NEON_HALF.vzero(), 1)),
                        assign(y, vload1_half(x[2].toPtr(NEON_HALF.t.t), u, 2))
                )))
	],
        # load using subvectors of length 2
        [(y,x,opts) -> assign(y, nth(x[1].toPtr(TVect(TReal, 2)), 0)) ]
    ],

    svstore := [
        [  (y,x,opts) -> vstore1_half(y[1].toPtr(NEON_HALF.t.t), x, 1),
           (y,x,opts) -> chain(
               vstore1_half(y[1].toPtr(NEON_HALF.t.t), x, 1),
               vstore1_half(y[2].toPtr(NEON_HALF.t.t), x, 2)
           )
        ],
        [  (y,x,opts) -> assign(nth(y[1].toPtr(NEON_HALF.t), 0), x) ]
    ],

    mul_cx := (self, opts) >> ((y,x,c) -> let(
	u := self.freshU(),  v1 := self.freshT(), v2 := self.freshT(), r1 := self.freshT(),
	decl([u, v1, v2, r1], chain(
		assign(r1, self.t.value([-1.0,1.0])),
         	assign(u, vtrnq_half(c, c)),
		assign(v1, x * vextract_half(u, [0])),
		assign(v2, vswapcx_half(x) * vextract_half(u, [1])),
		assign(y, fma(v1, v2, r1))))
    )),

    mul_cx_conj := (self, opts) >> ((y,x,c) -> let(
	u := self.freshU(),  v1 := self.freshT(), v2 := self.freshT(), r1 := self.freshT(),
	decl([u, v1, v2, r1], chain(
		assign(r1, self.t.value([1.0, -1.0])),
         	assign(u, vtrnq_half(c, c)),
		assign(v1, x * vextract_half(u, [0])),
		assign(v2, vswapcx_half(x) * vextract_half(u, [1])),
		assign(y, fma(v1, v2, r1))))
    )),

    swap_cx := (y, x, opts) -> assign(y, vswapcx_half(x)),
));

Class(NEON, SIMD_NEON, rec(
    info := "NEON 4 x floating point",

    countrec := rec(
        ops := [
            [ add, sub ],
	    [ mul, vmulcx_neon ],
            [ fma, fms, nfma ],
            [ vpacklo_neon,      vpackhi_neon,       vunpacklo_neon,    vunpackhi_neon,
              vunpacklolo2_neon, vunpacklohi2_neon,  vunpackhilo2_neon, vunpackhihi2_neon,
              vtransposelo_neon, vtransposehi_neon,  vrev_neon,         vswapcx_neon,
              vuzpq_32f,         vzipq_32f,          vtrnq_32f,         vextract_neon_4x32f,
	      vdup_lane_neon ],
            [ vload1_neon, vload_half_neon, vcombine_neon ],
            [ deref, nth, ],
            Value      # Value without [] is a keyword in countOps !!
        ],
        printstrings := ["[adds]","[mults]","[fmas]","[vperms]","[svldst]","[vldst]","[vval]"],
        arithcost := (self, opcount) >> opcount[1]+opcount[2]
    ),

    v     := 4,
    t     := TVect(T_Real(32), 4),
    freshU := self >> var.fresh_t("u", TVect(self.t, 2)),
    freshT := self >> var.fresh_t("t", self.t),

    instr := [vpacklo_neon,    vpackhi_neon,
              vunpacklo_neon,  vunpackhi_neon,
              vunpacklolo2_neon,   vunpacklohi2_neon,
              vunpackhilo2_neon,   vunpackhihi2_neon,
              vtransposelo_neon, vtransposehi_neon,
	      vrev_neon ],

    includes := self >> self.commonIncludes() ::
        ["<arm_neon.h>", "<include/omega32_neon.h>", "<include/mm_malloc.h>"],

    dupload  := (y, x) -> assign(y, vdup(x, 4)),
    duploadn := (y, x, n) -> assign(y, vdup_lane_neon(x, n)),

    reverse := (y,x) -> assign(vref(y,0,4), vrev_neon(vext_neon(vref(x,0,4), vref(x,0,4), 2))),
    RCVIxJ2 := (y,x,opts) -> assign(y, vext_neon(x, x, 2)),

    # this is shift "right" in Spiral notation (and ARM too)
    # [a b c d] -> [0 a b c]
    bin_shl1 := (self, y,x,opts) >> assign(y, vext_neon(self.t.zero(), x, 1)),
    # this is shift "left" in Spiral notation (and ARM too)
    # [a b c d] -> [b c d 0]
    bin_shr1 := (self, y,x,opts) >> assign(y, vext_neon(x, self.t.zero(), 3)),

    # this is shift "right" in Spiral notation (and ARM too)
    # [a b c d], [e f g h] -> [d e f g]
    bin_shl2 := (self, y,x,opts) >> assign(y, vext_neon(x[1], x[2], 1)),
    # this is shift "left" in Spiral notation (and ARM too)
    # [a b c d], [e f g h] -> [b c d e]
    bin_shr2 := (self, y,x,opts) >> assign(y, vext_neon(x[1], x[2], 3)),

    # [a b c d] [e f g h] -> [e d c b]
    # support for VO1dsJ(n, v)
    bin_shrev := (self, y, x, opts) >> let(
	u := var.fresh_t("T", self.t),
	decl([u], chain(
		assign(u, vext_neon(x[1], x[2], 3)), # [ b c d e]
		assign(u, vext_neon(u, u, 2)),       # [ d e b c]
		assign(y, vrev_neon(u))))            # [ e d c b]
    ),

    loadCont := (self, n, y, yofs, x, xofs, xofs_align, opts) >> let(
	nn := _unwrap(n),
	yy := vtref(self.t, y, yofs),
	p0 := nth(x, xofs).toPtr(self.t.t),
	p2 := nth(x, xofs+2).toPtr(self.t.t),
	Cond(
	    nn=1,
	        assign(yy, vload1_neon(p0, self.vzero(), 1)),
	    nn=2,
                assign(yy, vcombine_neon(
			       vload_half_neon(p0),
			       NEON_HALF.vzero())),
	    nn=3,
                assign(yy, vcombine_neon(
			       vload_half_neon(p0),
			       vload1_half(p2, NEON_HALF.vzero(), 1))),
	    nn=4,
	        assign(yy, deref(nth(x, xofs).toPtr(self.t)))
	)
    ),

    storeCont := (self, n, y, yofs, yofs_align, x, xofs, opts) >> let(
	nn := _unwrap(n),
	xx := vtref(self.t, x, xofs),
	p0 := nth(y, yofs).toPtr(self.t.t),
	p2 := nth(y, yofs+2).toPtr(self.t.t),
	Cond(
	    nn=1,
		vstore1_neon(p0, xx, 1),
	    nn=2,
		vstore2lo_neon(p0, xx),
	    nn=3, chain(
		vstore2lo_neon(p0, xx),
		vstore1_neon(p2, xx, 3)),
	    nn=4,
	        assign(deref(nth(y, yofs).toPtr(self.t)), xx)
	)
    ),

    svload := [
        # load using subvectors of length 1
	[
            (y,x,opts) -> assign(y, vload1_neon(x[1].toPtr(NEON.t.t), NEON.vzero(), 1)),

            (y,x,opts) -> let(u := NEON.freshT(),
                    decl([u], chain(
                            assign(u, vload1_neon(x[1].toPtr(NEON.t.t), NEON.vzero(), 1)),
                            assign(y, vload1_neon(x[2].toPtr(NEON.t.t), u,            2))))),
            (y,x,opts) -> let(u := NEON.freshT(),
                    decl([u], chain(
                            assign(u, vload1_neon(x[1].toPtr(NEON.t.t), NEON.vzero(), 1)),
                            assign(u, vload1_neon(x[2].toPtr(NEON.t.t), u,            2)),
                            assign(y, vload1_neon(x[3].toPtr(NEON.t.t), u,            3))))),
            (y,x,opts) -> let(u := NEON.freshT(),
                    decl([u], chain(
                            assign(u, vload1_neon(x[1].toPtr(NEON.t.t), NEON.vzero(), 1)),
                            assign(u, vload1_neon(x[2].toPtr(NEON.t.t), u,            2)),
                            assign(u, vload1_neon(x[3].toPtr(NEON.t.t), u,            3)),
                            assign(y, vload1_neon(x[4].toPtr(NEON.t.t), u,            4)))))
	],
        # load using subvectors of length 2
        [
            (y,x,opts) -> let(
		u := var.fresh_t("T", NEON_HALF.t),
                decl([u], chain(
                        assign(u, vload_half_neon(x[1].toPtr(NEON_HALF.t.t))),
                        assign(y, vcombine_neon(u, NEON_HALF.vzero()))))),
            (y,x,opts) -> let(
		u1 := var.fresh_t("T", NEON_HALF.t),
		u2 := var.fresh_t("T", NEON_HALF.t),
                decl([u1, u2] , chain(
                        assign(u1, vload_half_neon(x[1].toPtr(NEON_HALF.t.t))),
                        assign(u2, vload_half_neon(x[2].toPtr(NEON_HALF.t.t))),
			assign(y, vcombine_neon(u1, u2)))))
        ]],

    svstore := [
        # store using subvectors of length 1
        [  (y,x,opts) -> vstore1_neon(y[1].toPtr(NEON.t.t), x, 1),
           (y,x,opts) -> chain(
               vstore1_neon(y[1].toPtr(NEON.t.t), x, 1),
               vstore1_neon(y[2].toPtr(NEON.t.t), x, 2)),
           (y,x,opts) -> chain(
               vstore1_neon(y[1].toPtr(NEON.t.t), x, 1),
               vstore1_neon(y[2].toPtr(NEON.t.t), x, 2),
               vstore1_neon(y[3].toPtr(NEON.t.t), x, 3)),
           (y,x,opts) -> chain(
               vstore1_neon(y[1].toPtr(NEON.t.t), x, 1),
               vstore1_neon(y[2].toPtr(NEON.t.t), x, 2),
               vstore1_neon(y[3].toPtr(NEON.t.t), x, 3),
               vstore1_neon(y[4].toPtr(NEON.t.t), x, 4))
        ],
        # store using subvectors of length 2
        [  (y,x,opts) -> vstore2lo_neon(y[1].toPtr(NEON.t.t), x),
           (y,x,opts) -> chain(vstore2lo_neon(y[1].toPtr(NEON.t.t), x),
                               vstore2hi_neon(y[2].toPtr(NEON.t.t), x))
        ]
    ],

    mul_cx := (self, opts) >> ((y,x,c) -> let(
	    u := self.freshU(),  v1 := self.freshT(), v2 := self.freshT(), r1 := self.freshT(),
	    decl([u, v1, v2, r1], chain(
		    assign(r1, self.t.value([-1.0,1.0,-1.0,1.0])),
         	    assign(u, vtrnq_32f(c, c)),
		    assign(v1, x * vextract_neon_4x32f(u, [0])),
		    assign(v2, vswapcx_neon(x) * vextract_neon_4x32f(u, [1])),
		    assign(y, fma(v1, v2, r1))))
    )),

    mul_cx_conj := (self, opts) >> ((y,x,c) -> let(
	    u := self.freshU(),  v1 := self.freshT(), v2 := self.freshT(), r1 := self.freshT(),
	    decl([u, v1, v2, r1], chain(
		    assign(r1, self.t.value([1.0, -1.0, 1.0, -1.0])),
         	    assign(u, vtrnq_32f(c, c)),
		    assign(v1, x * vextract_neon_4x32f(u, [0])),
		    assign(v2, vswapcx_neon(x) * vextract_neon_4x32f(u, [1])),
		    assign(y, fma(v1, v2, r1))))
    )),

    swap_cx := (y, x, opts) -> assign(y, vswapcx_neon(x)),
));

SIMD_ISA_DB.addISA(NEON);
SIMD_ISA_DB.addISA(NEON_HALF);


text:
#
# Read("~/Workspace/groupsSB/epi/B2/handle.gi");
#

type:="B";
rank:=2;
nr_pos_roots:=4;

Read("~/Workspace/groupsSB/epi/group.gi");

#sList(rels,r->Value(r,[vars[22]],[vars[2]+vars[12]]);

#
# c_1, c_2 not 0 => c_3=c_1/2
#
# a test that the coeffs are right
#a:=root_group(1,xvars[1]);
#b:=root_group(2,xvars[2]);
#c:=root_group(3,xvars[3]);
#d:=root_group(4,xvars[4]);
#test1:=evaluate_U(a*b*c*d,[[[xvars[1]],[One(APR)]],[[xvars[2]],[One(APR)]],[[xvars[3]],[-One(APR)/2]],[[xvars[4]],[-2*One(APR)/3]]]);
#test2:=evaluate_U(a*b*c*d,[[[xvars[1]],[One(APR)*2]],[[xvars[2]],[One(APR)*2]],[[xvars[3]],[-One(APR)*2]],[[xvars[4]],[-16*One(APR)/3]]]);
#Set(Concatenation(test2-test1^2));
handleReg:=function()
	local ua,ub,uab,uabb,u,e,nn,vals,rels,i,v;
	ua:=root_group(1,1);
	ub:=root_group(2,1);
	uab:=root_group(3,-1/2);
	uabb:=root_group(4,-2/3);
	u:=ua*ub*uab*uabb;
	e:=u-id_mat;

	nn:=List([1..2*nr_pos_roots+rank],i->xvars[i]);

	vals:=[];
 	Append(vals,[[[xvars[6]],[Zero(APR)]]]);
# 	Append(vals,[[[xvars[4]],[Zero(APR)]]]);
  	Append(vals,[[[xvars[5]],[Zero(APR)]]]);
  	Append(vals,[[[xvars[7]],[Zero(APR)]]]);
  	Append(vals,[[[xvars[8]],[Zero(APR)]]]);
  	Append(vals,[[[xvars[9]],[2*xvars[10]]]]);
 	Append(vals,[[[xvars[10]],[Zero(APR)]]]);
  	Append(vals,[[[xvars[2]],[xvars[1]]]]);
  	Append(vals,[[[xvars[3]],[Zero(APR)]]]);
#  	Append(vals,[[[xvars[3]],[-2*xvars[1]/3]]]);
#  	##
	#Append(vals,[[[xvars[1]],[Zero(APR)]]]);
	#Append(vals,[[[xvars[3]],[Zero(APR)]]]);
	#Append(vals,[[[xvars[7]],[Zero(APR)]]]);

	#
	# atentie la ordinea transpunerilor
	#
	rels:=nn*TransposedMat(u)-nn;
	rels:=evaluate_rels(rels,vals);
	nn:=evaluate_rels(nn,vals);

	return [u,nn,rels,vals];
end;

#
# Done
#
handleRegModule:=function()
	local ua,ub,uab,uabb,u,e,nn,vals,rels,i,v;
	ua:=root_group(1,1);
	ub:=root_group(2,1);
	uab:=root_group(3,-1/2);
	uabb:=root_group(4,-2/3);
	u:=ua*ub*uab*uabb;

	nn:=One(APR)*[0,0,0,0,xvars[1],xvars[2],xvars[3],0,0,0];
	#
	# atentie la ordinea transpunerilor
	#
	nn:=nn*TransposedMat(u);

	rels:=[];
	Append(rels,[nn[1]-nn[2]]);
	# Append(rels,[nn[3]-2*nn[1]/3]);

	vals:=[];
	Append(vals,[[[xvars[3]],[-6*xvars[1]/5+8*xvars[2]/5]]]); # because I want the projection on U3 to be trivial
	Append(vals,[[[xvars[2]],[3*xvars[1]/4]]]); # because I want the projections on U1 and U2 to be equal
#	Append(vals,[[[xvars[3]],[13*xvars[1]/2-2*xvars[2]]]]);
#	Append(vals,[[[xvars[2]],[203*3*xvars[1]/(4*85)]]]);

	rels:=evaluate_rels(rels,vals);
	nn:=evaluate_rels(nn,vals);

	return [u,nn,rels,vals];
end;

#
# c_2=0 c_3=0
#
#
#
handle23:=function()
	local ua,ub,uab,uabb,u,e,nn,vals,rels,i,v;
	ua:=root_group(1,1);
	uabb:=root_group(4,1);
	u:=ua*uabb;
	id_mat:=u^0;
	e:=u-id_mat;

	nn:=List([1..2*nr_pos_roots+rank],i->xvars[i]);

	vals:=[];
	Append(vals,[[[xvars[7]],[Zero(APR)]]]);
	Append(vals,[[[xvars[8]],[Zero(APR)]]]);
	Append(vals,[[[xvars[5]],[Zero(APR)]]]);
	Append(vals,[[[xvars[10]],[Zero(APR)]]]); # for p>2
	Append(vals,[[[xvars[9]],[Zero(APR)]]]); # for p>2
	Append(vals,[[[xvars[2]],[-xvars[6]]]]);
	# ###
	# Append(vals,[[[xvars[1]],[Zero(APR)]]]);
	# Append(vals,[[[xvars[3]],[Zero(APR)]]]);
	# Append(vals,[[[xvars[4]],[Zero(APR)]]]);
	# Append(vals,[[[xvars[8]],[Zero(APR)]]]);
	# Append(vals,[[[xvars[9]],[Zero(APR)]]]);

	#
	# atentie la ordinea transpunerilor
	#
	rels:=nn*TransposedMat(u)-nn;
	rels:=evaluate_rels(rels,vals);
	nn:=evaluate_rels(nn,vals);

	return [u,nn,rels,vals];
end;

#
# c_2=0 c_4=0
#
handle24:=function()
	local ua,ub,uab,uabb,u,e,nn,vals,rels,i,v;
	ua:=root_group(1,1);
	uab:=root_group(3,1);
	u:=ua*uab;
	id_mat:=u^0;
	e:=u-id_mat;

	nn:=List([1..2*nr_pos_roots+rank],i->xvars[i]);

	vals:=[];
	Append(vals,[[[xvars[7]],[Zero(APR)]]]);
	Append(vals,[[[xvars[8]],[Zero(APR)]]]);
	Append(vals,[[[xvars[5]],[Zero(APR)]]]);
	Append(vals,[[[xvars[2]],[Zero(APR)]]]); # for p>2
	Append(vals,[[[xvars[9]],[Zero(APR)]]]);
	Append(vals,[[[xvars[10]],[-xvars[6]]]]); # for p>2
	# # ###
	# Append(vals,[[[xvars[1]],[Zero(APR)]]]);
	# Append(vals,[[[xvars[3]],[Zero(APR)]]]);
	# Append(vals,[[[xvars[4]],[Zero(APR)]]]);
	# Append(vals,[[[xvars[8]],[Zero(APR)]]]);
	# Append(vals,[[[xvars[9]],[Zero(APR)]]]);

	#
	# atentie la ordinea transpunerilor
	#
	rels:=nn*TransposedMat(u)-nn;
	rels:=evaluate_rels(rels,vals);
	nn:=evaluate_rels(nn,vals);

	return [u,nn,rels,vals];
end;

#
# c_2=0
#
handle2:=function()
	local ua,ub,uab,uabb,u,e,nn,vals,rels,i,v;
	ua:=root_group(1,1);
	uab:=root_group(3,1);
	uabb:=root_group(4,1);
	u:=ua*uab*uabb;
	id_mat:=u^0;
	e:=u-id_mat;

	nn:=List([1..2*nr_pos_roots+rank],i->xvars[i]);

	vals:=[];
	Append(vals,[[[xvars[5]],[xvars[8]]]]);
	Append(vals,[[[xvars[7]],[-xvars[8]]]]);
	Append(vals,[[[xvars[2]],[-xvars[6]+xvars[9]]]]);
	Append(vals,[[[xvars[6]],[xvars[9]-xvars[10]]]]); # only in char 2
	###
	# Append(vals,[[[xvars[1]],[Zero(APR)]]]);
	# Append(vals,[[[xvars[3]],[Zero(APR)]]]);
	# Append(vals,[[[xvars[4]],[Zero(APR)]]]);
	# Append(vals,[[[xvars[8]],[Zero(APR)]]]);
	# Append(vals,[[[xvars[9]],[Zero(APR)]]]);

	#
	# atentie la ordinea transpunerilor
	#
	rels:=nn*TransposedMat(u)-nn;
	for i in [1..Length(rels)] do
		for v in vals do
			Print(rels[i],"\n");
			rels[i]:=One(APR)*Value(rels[i],v[1],v[2]);
			nn[i]:=One(APR)*Value(nn[i],v[1],v[2]);
		od;
	od;

	return [u,nn,rels,vals];
end;

#
# regular
# c_1, c_2 not, c_3 not 0
#
handle1:=function()

end;

ext: gi

max_stars_repo_name: iuliansimion/groupsSB

max_issues_repo_name: iuliansimion/groupsSB

max_forks_repo_name: iuliansimion/groupsSB


text:
Y := function(f)
    local u;
    u := x -> x(x);
    return u(y -> f(a -> y(y)(a)));
end;

fib := function(f)
    local u;
    u := function(n)
        if n < 2 then
            return n;
        else
            return f(n-1) + f(n-2);
        fi;
    end;
    return u;
end;

Y(fib)(10);
# 55

fac := function(f)
    local u;
    u := function(n)
        if n < 2 then
            return 1;
        else
            return n*f(n-1);
        fi;
    end;
    return u;
end;

Y(fac)(8);
# 40320

ext: gap

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
#############################################################################
##
#W  automfam.gi              automgrp package                  Yevgen Muntyan
#W                                                             Dmytro Savchuk
##
#Y  Copyright (C) 2003 - 2018 Yevgen Muntyan, Dmytro Savchuk
##

###############################################################################
##
#R  IsAutomFamilyRep
##
##  Any object from category AutomFamily which is created here, lies in a
##  category IsAutomFamilyRep.
##  Family of Autom object <a> contains all essential information about
##  (mathematical) automaton which generates group containing <a>:
##  it contains automaton, properties of automaton and group generated by it,
##  etc.
##  Also family contains group generated by states of underlying automaton.
##
DeclareRepresentation("IsAutomFamilyRep",
                      IsComponentObjectRep and IsAttributeStoringRep,
                      [ "freegroup",      # IsAutom objects store words from this group
                        "freegens",       # list [f1, f2, ..., fn, f1^-1, ..., fn^-1, 1]
                                          # where f1..fn are generators of freegroup,
                                          # n is numstates; 1 is stored if and only if
                                          # trivstate is not zero.
                                          # Some fi^-1 may be missing if corresponding
                                          # generator is not invertible (but the list still
                                          # has the length of 2n+1).
                        "automgens",      # Generators of the group, list of length n.
                        "numstates",      # number of non-trivial generating states
                        "deg",
                        "trivstate",      # 0 or 2*numstates+1
                        "isgroup",        # whether all generators are invertible
                        "names",          # list of non-trivial generating states
                        "automatonlist",  # the automaton table, states correspond to freegens
                        "oldstates",      # mapping from states in the original table used to
                                          # define the group to the states in automatonlist:
                                          # Autom(fam!.freegens[oldtstates[k]], fam) is the element
                                          # which corresponds to k-th state in the original automaton
                        "rws",            # rewriting system
                        "use_rws",        # whether to use rewriting system in multiplication
                        "use_contraction" # whether to use contraction in IsOne
                      ]);

BindGlobal("AG_FixAutomList",
function(list, oldstates, names)
  local i, j, deg, isgroup, numstates, perm, trivstate;

  deg := Length(list[1]) - 1;
  isgroup := true;
  trivstate := 0;
  numstates := Length(list);

  for i in [1..Length(list)] do
    if AG_IsTrivialStateInList(i, list) then
      trivstate := i;
      numstates := numstates - 1;
      break;
    fi;
  od;

  # make sure trivial state in oldstates is right
  if trivstate <> 0 then
    for i in [1..Length(oldstates)] do
      if oldstates[i] = trivstate then
        oldstates[i] := 2*numstates + 1;
      elif oldstates[i] > trivstate then
        oldstates[i] := oldstates[i] - 1;
      fi;
    od;

    Remove(names, trivstate);
  fi;

  # move trivial state to the end
  if trivstate <> 0 then
    for i in [1..Length(list)] do
      for j in [1..deg] do
        if list[i][j] = trivstate then
          list[i][j] := 2*numstates + 1;
        elif list[i][j] > trivstate then
          list[i][j] := list[i][j] - 1;
        fi;
      od;
    od;

    list := Concatenation(list{[1..(trivstate-1)]},
                          list{[(trivstate+1)..Length(list)]});
    trivstate := 2*numstates + 1;
    list[trivstate] := List([1..deg], k -> trivstate);
    list[trivstate][deg+1] := ();
  fi;

  # add inverses of states
  for i in [1..numstates] do
    if AG_IsInvertibleStateInList(i, list) then
      list[i+numstates] := [];

      list[i][deg+1] := AG_PermFromTransformation(list[i][deg+1]);

      perm := list[i][deg+1];
      list[i+numstates][deg+1] := perm^-1;
      for j in [1..deg] do
        list[i+numstates][j] := list[i][j^(perm^-1)];
        if list[i+numstates][j] <> trivstate then
          list[i+numstates][j] := list[i+numstates][j] + numstates;
        fi;
      od;
    else
      isgroup := false;
      if list[i][deg+1]^-1<>fail then
        list[i][deg+1] := AG_PermFromTransformation(list[i][deg+1]);
      fi;
    fi;
  od;

  return [list, numstates, trivstate, isgroup];
end);

###############################################################################
##
#M  AutomFamily(<list>, <names>, <bind_vars>)
##
InstallMethod(AutomFamily, "for [IsList, IsList, IsBool]",
              [IsList, IsList, IsBool],
function (list, names, bind_global)
  local deg, tmp, trivstate, numstates, i,
        freegroup, freegens, a, family, oldstates,
        isgroup;

  if not AG_IsCorrectAutomatonList(list, false) then
    Error("in AutomFamily(IsList, IsList, IsString):\n",
          "  given list is not a correct list representing automaton\n");
  fi;

# 1. make a local copy of arguments, since they will be modified and put into the result

  list := StructuralCopy(list);
  names := StructuralCopy(names);
  deg := Length(list[1]) - 1;

# 2. Reduce automaton, find trivial state, permute states

  tmp := AG_ReducedAutomatonInList(list);
  list := tmp[1];
  oldstates := tmp[3];
  names := List(tmp[2], x -> names[x]);

  tmp := AG_FixAutomList(list, oldstates, names);
  list := tmp[1];
  numstates := tmp[2];
  trivstate := tmp[3];
  isgroup := tmp[4];

# 3. Create FreeGroup and FreeGens

  freegroup := FreeGroup(names);
  freegens := ShallowCopy(FreeGeneratorsOfFpGroup(freegroup));
  for i in [1..numstates] do
    if IsBound(list[i+numstates]) then
      freegens[i+numstates] := freegens[i]^-1;
    fi;
  od;
  if trivstate <> 0 then
    freegens[trivstate] := One(freegroup);
  fi;

# 4. Create family

  if isgroup then
    family := NewFamily("AutomFamily",
                        IsInvertibleAutom,
                        IsInvertibleAutom,
                        IsAutomFamily and IsAutomFamilyRep);
  else
    family := NewFamily("AutomFamily",
                        IsAutom,
                        IsAutom,
                        IsAutomFamily and IsAutomFamilyRep);
  fi;

  family!.isgroup := isgroup;
  family!.deg := deg;
  family!.numstates := numstates;
  family!.trivstate := trivstate;
  family!.names := names;
  family!.freegroup := freegroup;
  family!.freegens := freegens;
  family!.automatonlist := list;
  family!.oldstates := oldstates;
  family!.use_rws := false;
  family!.rws := fail;
  family!.use_contraction := false;

  SetIsActingOnBinaryTree(family, deg = 2);
  SetDegreeOfTree(family, deg);
  SetTopDegreeOfTree(family, deg);

  family!.automgens := [];
  for i in [1..Length(list)] do
    if IsBound(list[i]) then
      family!.automgens[i] :=
        __AG_CreateAutom(family, freegens[i],
                        List([1..deg], j -> freegens[list[i][j]]),
                        list[i][deg+1],
                        i > numstates or IsBound(list[i+numstates]));
    fi;
  od;

  # XXX It's evil to bind global names, consider AssignGeneratorVariables
  # XXX Check whether names are actually valid names for variables
  if bind_global then
    for i in [1..family!.numstates] do
      if IsBoundGlobal(family!.names[i]) then
        UnbindGlobal(family!.names[i]);
      fi;
      BindGlobal(family!.names[i], family!.automgens[i]);
      MakeReadWriteGlobal(family!.names[i]);
    od;
    #BindGlobal(AG_Globals.identity_symbol, One(family));
    #MakeReadWriteGlobal(AG_Globals.identity_symbol);
  fi;

  return family;
end);

###############################################################################
##
#M  AutomFamily(<list>)
##
InstallMethod(AutomFamily, "for [IsList]", [IsList],
function(list)
  return AutomFamily(list, false);
end);

###############################################################################
##
#M  AutomFamily(<list>, <bind_vars>)
##
InstallMethod(AutomFamily, "for [IsList, IsBool]", [IsList, IsBool],
function(list, bind_vars)
  if not AG_IsCorrectAutomatonList(list, false) then
    Error("in AutomFamily(IsList):\n",
          "  given list is not a correct list representing automaton\n");
  fi;

  return AutomFamily(list,
                     List([1..Length(list)],
                          i -> Concatenation(AG_Globals.state_symbol, String(i))),
                     bind_vars);
end);

###############################################################################
##
#M  AutomFamily(<list>, <names>)
##
InstallMethod(AutomFamily, [IsList, IsList],
function(list, names)
  if not AG_IsCorrectAutomatonList(list, false) then
    Error("in AutomFamily(IsList, IsList):\n",
          "  given list is not a correct list representing automaton\n");
  fi;

  return AutomFamily(list, names, AG_Globals.bind_vars_autom_family);
end);

###############################################################################
##
#M  DualAutomFamily(<fam>)
##
InstallMethod(DualAutomFamily, "for [IsAutomFamily]",
              [IsAutomFamily],
function(fam)
  local list, dual;
  list := [1..fam!.numstates];
  if fam!.trivstate <> 0 then
    Add(list, fam!.trivstate);
  fi;
  list := AG_MinimalSubAutomatonInlist(list, fam!.automatonlist)[1];
  if not AG_HasDualInList(list) then
    return 0;
  fi;
  dual := AutomFamily(AG_DualAutomatonList(list),
                      List([1..DegreeOfTree(fam)],
                           i -> Concatenation(AG_Globals.state_symbol_dual, String(i))));
  SetDualAutomFamily(dual, fam);
  return dual;
end);

###############################################################################
##
#M  One(<fam>)
##
InstallMethod(One, "for [IsAutomFamily]", [IsAutomFamily],
function(fam)
  return __AG_CreateAutom(fam, One(fam!.freegroup),
                         List([1..fam!.deg], i -> One(fam!.freegroup)),
                         (), true);
end);

###############################################################################
##
#M  GroupOfAutomFamily(<fam>)
##
InstallMethod(GroupOfAutomFamily, "for [IsAutomFamily]",
              [IsAutomFamily],
function(fam)
  local g;

  if not fam!.isgroup then
    return fail;
  fi;

  if fam!.numstates > 0 then
    g := GroupWithGenerators(fam!.automgens{[1..fam!.numstates]});
  else
    g := Group(One(fam));
  fi;

  SetUnderlyingAutomFamily(g, fam);
  SetIsGroupOfAutomFamily(g, true);
  # XXX
  SetGeneratingAutomatonList(g, fam!.automatonlist);
  SetAutomatonList(g, fam!.automatonlist);
  SetDegreeOfTree(g, fam!.deg);
  SetTopDegreeOfTree(g, fam!.deg);
  SetIsActingOnBinaryTree(g, fam!.deg = 2);
  SetIsAutomatonGroup(g, true);

  return g;
end);

###############################################################################
##
#M  SemigroupOfAutomFamily( <fam> )
##
InstallMethod(SemigroupOfAutomFamily, "for [IsAutomFamily]",
              [IsAutomFamily],
function(fam)
  local g;

  if fam!.trivstate <> 0 then
    if fam!.numstates = 0 then
      g := Group(One(fam!.automgens[fam!.trivstate]));
    else
      g := MonoidByGenerators(fam!.automgens{[1..fam!.numstates]});
    fi;
  else
    g := SemigroupByGenerators(fam!.automgens{[1..fam!.numstates]});
  fi;

  SetUnderlyingAutomFamily(g, fam);

  # XXX
  SetGeneratingAutomatonList(g, fam!.automatonlist);
  SetAutomatonList(g, fam!.automatonlist);

  SetDegreeOfTree(g, fam!.deg);
  SetTopDegreeOfTree(g, fam!.deg);
  SetIsActingOnBinaryTree(g, fam!.deg = 2);
  SetIsAutomatonSemigroup(g, true);

  return g;
end);

###############################################################################
##
#M  UnderlyingFreeMonoid( <G> )
##
InstallMethod(UnderlyingFreeMonoid, "for [IsAutomFamily]",
              [IsAutomFamily],
function(fam)
  local monoid;

  if fam!.numstates <> 0 then
    monoid := MonoidByGenerators(GeneratorsOfGroup(fam!.freegroup));
    SetSize(monoid, infinity);
  else
    monoid := MonoidByGenerators(fam!.freegens[1]);
    SetSize(monoid, 1);
  fi;

  return monoid;
end);

###############################################################################
##
#M  UnderlyingFreeGroup( <G> )
##
InstallMethod(UnderlyingFreeGroup, "for [IsAutomFamily]",
              [IsAutomFamily],
function(fam)
  return fam!.freegroup;
end);

###############################################################################
##
##  AG_AbelImagesGenerators(<fam>)
##
InstallMethod(AG_AbelImagesGenerators, "for [IsAutomFamily]",
              [IsAutomFamily],
function(fam)
  if not fam!.isgroup then
    Error("the group is not generated by an invertible automaton");
  fi;
  return AG_AbelImageAutomatonInList(fam!.automatonlist);
end);

###############################################################################
##
#M  DiagonalPowerOp(<fam>, <n>)
##
InstallMethod(DiagonalPowerOp, "for [IsAutomFamily, IsPosInt]",
              [IsAutomFamily, IsPosInt],
function(fam, n)
  local names, list, states, dlist;

  if not fam!.isgroup then
    Error("the group is not generated by an invertible automaton");
  fi;

  list := fam!.automatonlist;
  states := [1..fam!.numstates];
  names := fam!.names;
  if fam!.trivstate <> 0 then
    Add(states, fam!.trivstate);
    Add(names, AG_Globals!.identity_symbol);
  fi;
  dlist := AG_MinimalSubAutomatonInlist(states, list)[1];
  return AutomatonGroup(AG_DiagonalPowerInList(dlist, n, names)[1],
                    AG_DiagonalPowerInList(dlist, n, names)[2]);
end);

###############################################################################
##
#M  DiagonalPower(<obj>)
##
InstallOtherMethod(DiagonalPower, "for [IsObject]", [IsObject],
function(obj)
  return DiagonalPower(obj, 2);
end);

###############################################################################
##
#M  MultAutomAlphabet(<fam>, <n>)
##
InstallMethod(MultAutomAlphabetOp, "for [IsAutomFamily, IsPosInt]",
              [IsAutomFamily, IsPosInt],
function(fam, n)
  local list, states, dlist;
  list := fam!.automatonlist;
  states := [1..fam!.numstates];
  if fam!.trivstate <> 0 then Add(states, fam!.trivstate); fi;
  dlist := AG_MinimalSubAutomatonInlist(states, list)[1];
  return AutomatonGroup(AG_MultAlphabetInList(dlist, n));
end);

###############################################################################
##
#M  MultAutomAlphabet(<obj>)
##
InstallOtherMethod(MultAutomAlphabet, "for [IsObject]", [IsObject],
function(obj)
  return MultAutomAlphabet(obj, 2);
end);

#############################################################################
##
#M  GeneratorsOfOrderTwo(<fam>)
##
InstallMethod(GeneratorsOfOrderTwo, "for [IsAutomFamily]", [IsAutomFamily],
function(fam)
  if not fam!.isgroup then
    Error("the group is not generated by an invertible automaton");
  fi;
  return Filtered(GeneratorsOfGroup(GroupOfAutomFamily(fam)), g -> IsOne(g^2));
end);

#E

max_stars_repo_name: gap-packages/automgrp

max_issues_repo_name: gap-packages/automgrp

max_forks_repo_name: gap-packages/automgrp


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

ImportAll(paradigms.vector);

Class(LSKernel, SumsBase, BaseContainer, rec(
    doNotMarkBB:= true,
    abbrevs := [ ch -> [ch,0], (ch,ops) -> [ch,ops] ],
    new := (self, ch, ops) >> SPL(WithBases(self, rec(
        info := Cond(
            IsInt(ops) or IsRat(ops)
                or IsValue(ops) or IsExp(ops), rec(
                    opcount := When(IsValue(ops),ops.v,ops),
                    free := Set([]),
                    loadFunc := fId(ch.dimensions[2]),
                    storeFunc := fId(ch.dimensions[1])
                ),

            IsRec(ops), ops,

            Error("unknown info")
        ),

        _children := [ch],

        dimensions := ch.dimensions
    ))),

    rChildren := self >> [self._children[1], self.info],

    rSetChild := meth(self, n, what)
        if n=2 then self.info := what;
        elif n=1 then self._children[1] := what;
        else Error("<n> must be in [1..2]");
        fi;
    end,

    mergeInfo := (self, r1, r2) >> rec(
        opcount := r1.opcount + r2.opcount,
        free := Concat(r1.free, r2.free),
        loadFunc := r2.loadFunc,
        storeFunc := r1.storeFunc
    ),

    print := meth(self, indent, indentStep)
        local s,ch,first,newline;

        ch := [self.child(1),self.info];
        if self._short_print or ForAll(ch, x->IsSPLSym(x) or IsSPLMat(x)) then
            newline := Ignore;
        else
            newline := self._newline;
        fi;

        first:=true;
        Print(self.__name__, "(");
        for s in ch do
            if(first) then first:=false;
            else Print(", "); fi;
            newline(indent + indentStep);
            When(IsSPL(s) or (IsRec(s) and IsBound(s.print) and NumGenArgs(s.print)=2),
                s.print(indent + indentStep, indentStep), Print(s));
        od;
        newline(indent);
        Print(")");
        self.printA();
    end,
));

Class(DMAGath, Gath, rec(doNotMarkBB:= true));

Class(DMAScat, Scat, rec(doNotMarkBB:= true));

Class(DMAFence, Buf, rec(doNotMarkBB:= true));

Class(SWPSum, ISum, rec(doNotMarkBB := true));

Class(DMAGathV, VGath, rec(doNotMarkBB := true));
Class(DMAScatV, VScat, rec(doNotMarkBB := true));


text:
# Built-in
Combinations([1 .. n], m);

Combinations([1 .. 5], 3);
# [ [ 1, 2, 3 ], [ 1, 2, 4 ], [ 1, 2, 5 ], [ 1, 3, 4 ], [ 1, 3, 5 ],
#   [ 1, 4, 5 ], [ 2, 3, 4 ], [ 2, 3, 5 ], [ 2, 4, 5 ], [ 3, 4, 5 ] ]


text:
# Code must be in 0 .. 255.
CharInt(65);
# 'A'
IntChar('Z');
# 90


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

#   count number of vector instructions in a ruletree
#   used to build vector base cases with minimal instruction count

VCost := function(R, O)
   return SPLRuleTree(R).vcost();
end;

Tensor.vcost := self >> let(C := self.children(), When(IsIdentitySPL(C[1]), C[1].params*C[2].vcost(), 1E15));
VTensor.vcost := self >> 0;
Compose.vcost := self >> Sum(List(self.children(), i->i.vcost()));
BlockVPerm.vcost := self >> self.n * self.child(1).vcost();


text:
%Terminals
    u0000
    u0001
    u0002
    u0003
    u0004
    u0005
    u0006
    u0007
    u0008
    u0009 -- HT
    u000A -- LF
    u000B
    u000C -- FF
    u000D -- CR
    u000E
    u000F
    u0010
    u0011
    u0012
    u0013
    u0014
    u0015
    u0016
    u0017
    u0018
    u0019
    u001A
    u001B
    u001C
    u001D
    u001E
    u001F
    u0020 ::= ' '
    u0021 ::= '!'
    u0022 ::= '"'
    u0023 ::= '%'
    u0024 ::= '%'
    u0025 ::= '%'
    u0026 ::= '&'
    u0027 ::= "'"
    u0028 ::= '('
    u0029 ::= ')'
    u002A ::= '*'
    u002B ::= '+'
    u002C ::= ','
    u002D ::= '-'
    u002E ::= '.'
    u002F ::= '/'
    u0030 ::= '0'
    u0031 ::= '1'
    u0032 ::= '2'
    u0033 ::= '3'
    u0034 ::= '4'
    u0035 ::= '5'
    u0036 ::= '6'
    u0037 ::= '7'
    u0038 ::= '8'
    u0039 ::= '9'
    u003A ::= ':'
    u003B ::= ';'
    u003C ::= '<'
    u003D ::= '='
    u003E ::= '>'
    u003F ::= '?'
    u0040 ::= '@'
    u0041 ::= 'A'
    u0042 ::= 'B'
    u0043 ::= 'C'
    u0044 ::= 'D'
    u0045 ::= 'E'
    u0046 ::= 'F'
    u0047 ::= 'G'
    u0048 ::= 'H'
    u0049 ::= 'I'
    u004A ::= 'J'
    u004B ::= 'K'
    u004C ::= 'L'
    u004D ::= 'M'
    u004E ::= 'N'
    u004F ::= 'O'
    u0050 ::= 'P'
    u0051 ::= 'Q'
    u0052 ::= 'R'
    u0053 ::= 'S'
    u0054 ::= 'T'
    u0055 ::= 'U'
    u0056 ::= 'V'
    u0057 ::= 'W'
    u0058 ::= 'X'
    u0059 ::= 'Y'
    u005A ::= 'Z'
    u005B ::= '['
    u005C ::= '\'
    u005D ::= ']'
    u005E ::= '^'
    u005F ::= '_'
    u0060 ::= '`'
    u0061 ::= 'a'
    u0062 ::= 'b'
    u0063 ::= 'c'
    u0064 ::= 'd'
    u0065 ::= 'e'
    u0066 ::= 'f'
    u0067 ::= 'g'
    u0068 ::= 'h'
    u0069 ::= 'i'
    u006A ::= 'j'
    u006B ::= 'k'
    u006C ::= 'l'
    u006D ::= 'm'
    u006E ::= 'n'
    u006F ::= 'o'
    u0070 ::= 'p'
    u0071 ::= 'q'
    u0072 ::= 'r'
    u0073 ::= 's'
    u0074 ::= 't'
    u0075 ::= 'u'
    u0076 ::= 'v'
    u0077 ::= 'w'
    u0078 ::= 'x'
    u0079 ::= 'y'
    u007A ::= 'z'
    u007B ::= '{'
    u007C ::= '|'
    u007D ::= '}'
    u007E ::= '~'
    u007F

    UNUSED
%End
%Trailers
/.
        //
        //
        //
       export   class %super_stream_class extends LpgLexStream
        {

         //
        //
        //
        static   tokenKind : number[]=  new Array(0x10000)  ; // 0x10000 == 65536
        static  __b_init : boolean = %super_stream_class.init_block(%super_stream_class.tokenKind);
        static  init_block(tokenKind : number[]) : boolean
        {
            for (let i = 0; i < tokenKind.length; ++i) {
                tokenKind[i] = 0;
            }
            tokenKind[0x0000] = %sym_type.%prefix%u0000%suffix%;           // 000    0x00
            tokenKind[0x0001] = %sym_type.%prefix%u0001%suffix%;           // 001    0x01
            tokenKind[0x0002] = %sym_type.%prefix%u0002%suffix%;           // 002    0x02
            tokenKind[0x0003] = %sym_type.%prefix%u0003%suffix%;           // 003    0x03
            tokenKind[0x0004] = %sym_type.%prefix%u0004%suffix%;           // 004    0x04
            tokenKind[0x0005] = %sym_type.%prefix%u0005%suffix%;           // 005    0x05
            tokenKind[0x0006] = %sym_type.%prefix%u0006%suffix%;           // 006    0x06
            tokenKind[0x0007] = %sym_type.%prefix%u0007%suffix%;           // 007    0x07
            tokenKind[0x0008] = %sym_type.%prefix%u0008%suffix%;           // 008    0x08
            tokenKind[0x0009] = %sym_type.%prefix%HT%suffix%;              // 009    0x09
            tokenKind[0x000A] = %sym_type.%prefix%LF%suffix%;              // 010    0x0A
            tokenKind[0x000B] = %sym_type.%prefix%u000B%suffix%;           // 011    0x0B
            tokenKind[0x000C] = %sym_type.%prefix%FF%suffix%;              // 012    0x0C
            tokenKind[0x000D] = %sym_type.%prefix%CR%suffix%;              // 013    0x0D
            tokenKind[0x000E] = %sym_type.%prefix%u000E%suffix%;           // 014    0x0E
            tokenKind[0x000F] = %sym_type.%prefix%u000F%suffix%;           // 015    0x0F
            tokenKind[0x0010] = %sym_type.%prefix%u0010%suffix%;           // 016    0x10
            tokenKind[0x0011] = %sym_type.%prefix%u0011%suffix%;           // 017    0x11
            tokenKind[0x0012] = %sym_type.%prefix%u0012%suffix%;           // 018    0x12
            tokenKind[0x0013] = %sym_type.%prefix%u0013%suffix%;           // 019    0x13
            tokenKind[0x0014] = %sym_type.%prefix%u0014%suffix%;           // 020    0x14
            tokenKind[0x0015] = %sym_type.%prefix%u0015%suffix%;           // 021    0x15
            tokenKind[0x0016] = %sym_type.%prefix%u0016%suffix%;           // 022    0x16
            tokenKind[0x0017] = %sym_type.%prefix%u0017%suffix%;           // 023    0x17
            tokenKind[0x0018] = %sym_type.%prefix%u0018%suffix%;           // 024    0x18
            tokenKind[0x0019] = %sym_type.%prefix%u0019%suffix%;           // 025    0x19
            tokenKind[0x001A] = %sym_type.%prefix%u001A%suffix%;           // 026    0x1A
            tokenKind[0x001B] = %sym_type.%prefix%u001B%suffix%;           // 027    0x1B
            tokenKind[0x001C] = %sym_type.%prefix%u001C%suffix%;           // 028    0x1C
            tokenKind[0x001D] = %sym_type.%prefix%u001D%suffix%;           // 029    0x1D
            tokenKind[0x001E] = %sym_type.%prefix%u001E%suffix%;           // 030    0x1E
            tokenKind[0x001F] = %sym_type.%prefix%u001F%suffix%;           // 031    0x1F
            tokenKind[0x0020] = %sym_type.%prefix%u0020%suffix%;           // 032    0x20
            tokenKind[0x0021] = %sym_type.%prefix%u0021%suffix%;           // 033    0x21
            tokenKind[0x0022] = %sym_type.%prefix%u0022%suffix%;           // 034    0x22
            tokenKind[0x0023] = %sym_type.%prefix%u0023%suffix%;           // 035    0x23
            tokenKind[0x0024] = %sym_type.%prefix%u0024%suffix%;           // 036    0x24
            tokenKind[0x0025] = %sym_type.%prefix%u0025%suffix%;           // 037    0x25
            tokenKind[0x0026] = %sym_type.%prefix%u0026%suffix%;           // 038    0x26
            tokenKind[0x0027] = %sym_type.%prefix%u0027%suffix%;           // 039    0x27
            tokenKind[0x0028] = %sym_type.%prefix%u0028%suffix%;           // 040    0x28
            tokenKind[0x0029] = %sym_type.%prefix%u0029%suffix%;           // 041    0x29
            tokenKind[0x002A] = %sym_type.%prefix%u002A%suffix%;           // 042    0x2A
            tokenKind[0x002B] = %sym_type.%prefix%u002B%suffix%;           // 043    0x2B
            tokenKind[0x002C] = %sym_type.%prefix%u002C%suffix%;           // 044    0x2C
            tokenKind[0x002D] = %sym_type.%prefix%u002D%suffix%;           // 045    0x2D
            tokenKind[0x002E] = %sym_type.%prefix%u002E%suffix%;           // 046    0x2E
            tokenKind[0x002F] = %sym_type.%prefix%u002F%suffix%;           // 047    0x2F
            tokenKind[0x0030] = %sym_type.%prefix%u0030%suffix%;           // 048    0x30
            tokenKind[0x0031] = %sym_type.%prefix%u0031%suffix%;           // 049    0x31
            tokenKind[0x0032] = %sym_type.%prefix%u0032%suffix%;           // 050    0x32
            tokenKind[0x0033] = %sym_type.%prefix%u0033%suffix%;           // 051    0x33
            tokenKind[0x0034] = %sym_type.%prefix%u0034%suffix%;           // 052    0x34
            tokenKind[0x0035] = %sym_type.%prefix%u0035%suffix%;           // 053    0x35
            tokenKind[0x0036] = %sym_type.%prefix%u0036%suffix%;           // 054    0x36
            tokenKind[0x0037] = %sym_type.%prefix%u0037%suffix%;           // 055    0x37
            tokenKind[0x0038] = %sym_type.%prefix%u0038%suffix%;           // 056    0x38
            tokenKind[0x0039] = %sym_type.%prefix%u0039%suffix%;           // 057    0x39
            tokenKind[0x003A] = %sym_type.%prefix%u003A%suffix%;           // 058    0x3A
            tokenKind[0x003B] = %sym_type.%prefix%u003B%suffix%;           // 059    0x3B
            tokenKind[0x003C] = %sym_type.%prefix%u003C%suffix%;           // 060    0x3C
            tokenKind[0x003D] = %sym_type.%prefix%u003D%suffix%;           // 061    0x3D
            tokenKind[0x003E] = %sym_type.%prefix%u003E%suffix%;           // 062    0x3E
            tokenKind[0x003F] = %sym_type.%prefix%u003F%suffix%;           // 063    0x3F
            tokenKind[0x0040] = %sym_type.%prefix%u0040%suffix%;           // 064    0x40
            tokenKind[0x0041] = %sym_type.%prefix%u0041%suffix%;           // 065    0x41
            tokenKind[0x0042] = %sym_type.%prefix%u0042%suffix%;           // 066    0x42
            tokenKind[0x0043] = %sym_type.%prefix%u0043%suffix%;           // 067    0x43
            tokenKind[0x0044] = %sym_type.%prefix%u0044%suffix%;           // 068    0x44
            tokenKind[0x0045] = %sym_type.%prefix%u0045%suffix%;           // 069    0x45
            tokenKind[0x0046] = %sym_type.%prefix%u0046%suffix%;           // 070    0x46
            tokenKind[0x0047] = %sym_type.%prefix%u0047%suffix%;           // 071    0x47
            tokenKind[0x0048] = %sym_type.%prefix%u0048%suffix%;           // 072    0x48
            tokenKind[0x0049] = %sym_type.%prefix%u0049%suffix%;           // 073    0x49
            tokenKind[0x004A] = %sym_type.%prefix%u004A%suffix%;           // 074    0x4A
            tokenKind[0x004B] = %sym_type.%prefix%u004B%suffix%;           // 075    0x4B
            tokenKind[0x004C] = %sym_type.%prefix%u004C%suffix%;           // 076    0x4C
            tokenKind[0x004D] = %sym_type.%prefix%u004D%suffix%;           // 077    0x4D
            tokenKind[0x004E] = %sym_type.%prefix%u004E%suffix%;           // 078    0x4E
            tokenKind[0x004F] = %sym_type.%prefix%u004F%suffix%;           // 079    0x4F
            tokenKind[0x0050] = %sym_type.%prefix%u0050%suffix%;           // 080    0x50
            tokenKind[0x0051] = %sym_type.%prefix%u0051%suffix%;           // 081    0x51
            tokenKind[0x0052] = %sym_type.%prefix%u0052%suffix%;           // 082    0x52
            tokenKind[0x0053] = %sym_type.%prefix%u0053%suffix%;           // 083    0x53
            tokenKind[0x0054] = %sym_type.%prefix%u0054%suffix%;           // 084    0x54
            tokenKind[0x0055] = %sym_type.%prefix%u0055%suffix%;           // 085    0x55
            tokenKind[0x0056] = %sym_type.%prefix%u0056%suffix%;           // 086    0x56
            tokenKind[0x0057] = %sym_type.%prefix%u0057%suffix%;           // 087    0x57
            tokenKind[0x0058] = %sym_type.%prefix%u0058%suffix%;           // 088    0x58
            tokenKind[0x0059] = %sym_type.%prefix%u0059%suffix%;           // 089    0x59
            tokenKind[0x005A] = %sym_type.%prefix%u005A%suffix%;           // 090    0x5A
            tokenKind[0x005B] = %sym_type.%prefix%u005B%suffix%;           // 091    0x5B
            tokenKind[0x005C] = %sym_type.%prefix%u005C%suffix%;           // 092    0x5C
            tokenKind[0x005D] = %sym_type.%prefix%u005D%suffix%;           // 093    0x5D
            tokenKind[0x005E] = %sym_type.%prefix%u005E%suffix%;           // 094    0x5E
            tokenKind[0x005F] = %sym_type.%prefix%u005F%suffix%;           // 095    0x5F
            tokenKind[0x0060] = %sym_type.%prefix%u0060%suffix%;           // 096    0x60
            tokenKind[0x0061] = %sym_type.%prefix%u0061%suffix%;           // 097    0x61
            tokenKind[0x0062] = %sym_type.%prefix%u0062%suffix%;           // 098    0x62
            tokenKind[0x0063] = %sym_type.%prefix%u0063%suffix%;           // 099    0x63
            tokenKind[0x0064] = %sym_type.%prefix%u0064%suffix%;           // 100    0x64
            tokenKind[0x0065] = %sym_type.%prefix%u0065%suffix%;           // 101    0x65
            tokenKind[0x0066] = %sym_type.%prefix%u0066%suffix%;           // 102    0x66
            tokenKind[0x0067] = %sym_type.%prefix%u0067%suffix%;           // 103    0x67
            tokenKind[0x0068] = %sym_type.%prefix%u0068%suffix%;           // 104    0x68
            tokenKind[0x0069] = %sym_type.%prefix%u0069%suffix%;           // 105    0x69
            tokenKind[0x006A] = %sym_type.%prefix%u006A%suffix%;           // 106    0x6A
            tokenKind[0x006B] = %sym_type.%prefix%u006B%suffix%;           // 107    0x6B
            tokenKind[0x006C] = %sym_type.%prefix%u006C%suffix%;           // 108    0x6C
            tokenKind[0x006D] = %sym_type.%prefix%u006D%suffix%;           // 109    0x6D
            tokenKind[0x006E] = %sym_type.%prefix%u006E%suffix%;           // 110    0x6E
            tokenKind[0x006F] = %sym_type.%prefix%u006F%suffix%;           // 111    0x6F
            tokenKind[0x0070] = %sym_type.%prefix%u0070%suffix%;           // 112    0x70
            tokenKind[0x0071] = %sym_type.%prefix%u0071%suffix%;           // 113    0x71
            tokenKind[0x0072] = %sym_type.%prefix%u0072%suffix%;           // 114    0x72
            tokenKind[0x0073] = %sym_type.%prefix%u0073%suffix%;           // 115    0x73
            tokenKind[0x0074] = %sym_type.%prefix%u0074%suffix%;           // 116    0x74
            tokenKind[0x0075] = %sym_type.%prefix%u0075%suffix%;           // 117    0x75
            tokenKind[0x0076] = %sym_type.%prefix%u0076%suffix%;           // 118    0x76
            tokenKind[0x0077] = %sym_type.%prefix%u0077%suffix%;           // 119    0x77
            tokenKind[0x0078] = %sym_type.%prefix%u0078%suffix%;           // 120    0x78
            tokenKind[0x0079] = %sym_type.%prefix%u0079%suffix%;           // 121    0x79
            tokenKind[0x007A] = %sym_type.%prefix%u007A%suffix%;           // 122    0x7A
            tokenKind[0x007B] = %sym_type.%prefix%u007B%suffix%;           // 123    0x7B
            tokenKind[0x007C] = %sym_type.%prefix%u007C%suffix%;           // 124    0x7C
            tokenKind[0x007D] = %sym_type.%prefix%u007D%suffix%;           // 125    0x7D
            tokenKind[0x007E] = %sym_type.%prefix%u007E%suffix%;           // 126    0x7E
            tokenKind[0x007F] = %sym_type.%prefix%u007F%suffix%;           // 127    0x7F

            tokenKind[0xFFFF] = %sym_type.%prefix%EOF%suffix%;

            //
            // Every other character not yet assigned is treated initially as unused
            //
            for (let i = 0x007F; i < 0xFFFF; i++)
                if (tokenKind[i] == 0) tokenKind[i] = %sym_type.%prefix%UNUSED%suffix%;
            return true;
        }

        public    getKind(number i) : number // Classify character at ith location
        {
            return (i >= this.getStreamLength()
                       ? 0xffff
                       : %super_stream_class.tokenKind[getIntValue(i)]);
        }

        public  orderedExportedSymbols() : string[] { return %exp_type.orderedTerminalSymbols; }

        constructor(fileName: string, inputChars?: string, tab: number=4) {
             super(fileName, inputChars, tab);
         }
        }
./
%End
%Headers
    --
    -- Additional methods for the action class not provided in the template
    --
    /.
        //
        // The Lexer contains an array of characters as the input stream to be parsed.
        // There are methods to retrieve and classify characters.
        // The lexparser "token" is implemented simply as the index of the next character in the array.
        // The Lexer : the abstract class LpgLexStream with an implementation of the abstract
        // method getKind.  The template defines the Lexer class and the lexer() method.
        // A driver creates the action class, "Lexer", passing an Option object to the constructor.
        //
       kwLexer :  %kw_lexer_class;
       public   printTokens : boolean =false;

       private static  readonly   ECLIPSE_TAB_VALUE: number = 4;

        public  getKeywordKinds() : number[] { return this.kwLexer.getKeywordKinds(); }

        constructor(filename : string) : this(filename, ECLIPSE_TAB_VALUE)
        {

            this.kwLexer = new %kw_lexer_class(this.getInputChars(), %_IDENTIFIER);
        }

        public  initialize(filename : string,content? : string) : void
        {
            super.initialize(filename,content);
            if (this.kwLexer == null)
                 this.kwLexer = new %kw_lexer_class(this.getInputChars(), %_IDENTIFIER);
            else this.kwLexer.setInputChars(this.getInputChars());
        }

         void makeToken(kind : number)
        {
            let startOffset = this.getLeftSpan(),
                endOffset = this.getRightSpan();
            this.makeToken(startOffset, endOffset, kind);
            if (this.printTokens) this.printValue(startOffset, endOffset);
        }

         void makeComment(kind : number)
        {
            let startOffset = this.getLeftSpan(),
                endOffset = this.getRightSpan();
            super.getPrsStream().makeAdjunct(startOffset, endOffset, kind);
        }

         skipToken() : void
        {
            if (this.printTokens) this.printValue(this.getLeftSpan(), this.getRightSpan());
        }

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

         printValue(startOffset : number, endOffset : number) : void
        {
             let s = lexStream.getInputChars().substr(startOffset, endOffset - startOffset + 1);
             console.Out.Write(s);
        }

    ./
%End


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

Declare(param);

#Special ordering for parameters
#Allows control on the interfaces of autolib
Class(ParamOps, ExpOps, rec(
   \<   := (e1,e2) -> Cond(
       IsRec(e1) and IsRec(e2)
       and IsBound(e1.sort_weight)
       and IsBound(e2.sort_weight)
       and e1.sort_weight<>e2.sort_weight,
       e1.sort_weight < e2.sort_weight,
       ObjId(e1) <> ObjId(e2), ObjId(e1) < ObjId(e2),
       e1.rChildren() < e2.rChildren())
));

# Autolib coldness constants.
# pcCold -  parameter must be available at plan time;
# pcHot - parameter available at compute time only;
# pcAny - parameter can be hot or cold (default);
#

Class(pcCold, ConstClass);
Class(pcHot,  ConstClass);
Class(pcAny,  ConstClass);

#F param(<type>, <id>, [sort weight, [coldness]])  - symbolic representation of a parameter
#F                        (eg. of a codelet)
#F Fields: .t  - types
#F         .id - name of the parameter
#F         .sort_weight - weight during sorting
#F         .coldness - parameter coldness in autolib interfaces
Class(param, Loc, rec(
    isParam := true,
    __call__ := (arg) >> let(
                    self := arg[1],
                    type := arg[2],
                    id   := arg[3],
                    WithBases(self, Checked(IsString(id), IsType(type), rec(
                            operations := ParamOps,
                            id := id,
                            t := type,
                            sort_weight := When(Length(arg)>=4, arg[4], 0),
                            coldness := When(Length(arg)>=5, arg[5], pcAny))))),

    print := self >> Print(self.name, "(", self.t, ", \"", self.id, "\")"),
    rChildren := self >> [self.t, self.id],
    rSetChild := rSetChildFields("t", "id"),
    from_rChildren := (self, rch) >> CopyFields(ObjId(self)(rch[1], rch[2]),
            rec(sort_weight := self.sort_weight, coldness := self.coldness)),
    eval := self >> self,
    can_fold := False,

    computeType := self >> self.t,
    setRange := meth(self, r) # when used as a loop counter variable
       self.range := r;
       return self;
    end,
));

IsParam := x -> IsRec(x) and IsBound(x.isParam) and x.isParam;

Class(in_param, param);

Class(Unk, Exp, rec(
    __call__ := (self, t) >> Checked(IsType(t),
	WithBases(self, rec(operations := ExpOps, args:=[t], t := t))),
    computeType := self >> self.args[1],

   isUnk := true,
));

IsUnk := x -> IsRec(x) and IsBound(x.isUnk) and x.isUnk;

Class(UnkInt, Unk(TInt), rec(print := self >> Print(self.__name__)));

#F allocate(<var>, <type>) -- equivalent of malloc()
#F See also: deallocate(), zallocate()
Class(allocate, assign);

#F zero-allocate(<var>, <type>) -- equivalent of calloc()
#F See also: deallocate(), allocate()
Class(zallocate, assign);

#F deallocate(<var>, <type>) -- equivalent of free()
#F See also: allocate(), zallocate()
Class(deallocate, assign);

# fld(<type>, <loc>, <field_id>)
Class(fld, Loc, rec(
    __call__ := (self, type, loc, id) >> Checked(IsString(id), IsType(type), IsLoc(loc),
	WithBases(self, rec(
		operations := ExpOps,
		loc := loc,
		id := id,
		t := type))),
    computeType := self >> self.t,
    print := self >> Print(self.name, "(", self.t, ", ", self.loc, ", \"", self.id, "\")"),
    rChildren := self >> [self.t, self.loc, self.id],
    rSetChild := rSetChildFields("t", "loc", "id"),
    eval := self >> self,
    can_fold := False,
));

# ufld(<loc>, <field_id>)
Class(ufld, fld, rec(
    __call__ := (self, loc, id) >> Checked(IsString(id), IsLoc(loc),
	WithBases(self, rec(
		operations := ExpOps,
		loc := loc,
		id := id,
		t := TUnknown))),
    print := self >> Print(self.name, "(", self.loc, ", \"", self.id, "\")"),
    rChildren := self >> [self.loc, self.id],
    rSetChild := rSetChildFields("loc", "id"),
    eval := self >> self
));

# struct(<id>, <fields>)
#   <id> is the string that gives the name of the structure
#   <fields> is a list of <param>s or <var>s that go into the structure
#
Class(struct, Command, rec(
    __call__ := (self, id, fields) >> Checked(IsString(id), IsList(fields), WithBases(self,
	rec(operations := CmdOps,
	    id := id,
	    fields := fields))),
    rChildren := self >> [self.id, self.fields],
    rSetChild := rSetChildFields("id", "fields"),
    print := (self, i, si) >> Print(self.name, "(\"", self.id, "\", ", self.fields, ")"),
));

Class(ret, ExpCommand);

SubstParams := (s,bindings) ->
    SubstTopDownNR(s, @(1, param, e -> IsBound(bindings.(e.id))), e -> bindings.(e.id));

SubstParamsCustom := (s,bindings,objid_list) ->
    SubstTopDownNR(s, @(1, objid_list, e -> IsBound(bindings.(e.id))), e -> bindings.(e.id));

SubstVarsSafe := (s,bindings) ->
    SubstTopDownNR(s, @(1, var, e -> IsBound(bindings.(e.id))), e -> bindings.(e.id));


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

## The hoister is a compiler pass that separates looped code in three
## parts:
##  - the hoisted code a.k.a before the loop
##  - the body of the loop, a.k.a the loop itself
##  - the epilogue of the loop, a.k.a what is run after the loop
##
## Essentially, it is done by separating statements that depend on the loop
## index, that have to be in the body from statements that do not which can
## be hoisted.
##
## On top of this, the following code performs two optimizations:
##  - if code is being stored to a memory location that is loop independent,
##    then this memory location is scalarized, i.e. the store is performed in
##    register and the memory write is delayed (put inside the epilogue)
##  - The code performs loop variable induction. This replaces statements that
##    depend of param*i, by accumulators that add plus param at each pass of
##    the loop. The code also attempts to discover the minimal accumulation
##    value to have a few number of window pointers

Hoister := function(c)
    local loop_idx, loop_range, idxpool, hoist, body, epilogue, lvlist, lvfifo, asns, asn, lv, freevars,
    candidate_memory_acc, confirmed_memory_acc, stmt, stmts, lv, freevars, equal_second_member, d, step,
    lvfamilies, fams, fam, s, cc, pointer, j, var_usage, members, member_usage, rootexpr, m, newloop_idx, newloop_range, newloop_start, newloop_inc ;
    #Hoister only applies to loopn(decl(chain(...)))
    if ((ObjId(c)=loopn) and (ObjId(c.cmd)=decl) and (ObjId(c.cmd.cmd)=chain)) then

         c := Compile.pullDataDeclsRefs(c);
         loop_idx := c.var;
         loop_range := c.range;
         c := c.cmd;
         MarkDefUse(c);
         #idxpool will contain the set of all variables that change with the loop index
         #(and therefore should not be hoisted)
         idxpool := Set([loop_idx]);

         #Loop variable induction
         #This gathers all expressions that are linear in the index variable
         #and stores them as (v, start, inc) in the loop-variable list (lvlist)
         lvlist := [];
         lvfifo := [];

         Add(lvfifo, rec(v := loop_idx,
                     start := loop_idx.t.zero(),
                     inc := V(1)));

         while(Length(lvfifo)>0) do

             lv := lvfifo[1];
             lvfifo := ListWithout(lvfifo, 1);
             Add(lvlist, lv);

             #If lv.start is a variable, add this variable to the lvfifo
             if IsVar(lv.start) then
                  #If lv.start is used twice, it cannot be used as an induction variable.
                  if (Length(Collect(c, [assign, lv.start, @(0)]))=1
                      and Length(Collect(c, lv.start))=1) then
                      SubstBottomUp(c, lv.v, e->lv.start);
                      SubstBottomUp(lvfifo, lv.v, e->lv.start);
                      SubstBottomUp(lvlist, lv.v, e->lv.start);

                      [asns, c] := Pull(c, [assign, lv.start, @(0)],e->skip(), e->e);

                      Add(lvfifo, rec(v := lv.start,
                              start := asns[1].exp,
                              inc := lv.inc));
                      AddSet(idxpool, lv.start);

                 fi;
             fi;

             #Gather all mults and adds that depend on lv and add them too
             [asns, c] := Pull(c, [assign, @(0), [add, @(1).cond(e->ObjId(e)<>Value), lv.v]],
                 e->skip(), e->e);
             for asn in asns do
                 Add(lvfifo, rec(v := asn.loc,
                         start := add(asn.exp.args[1], lv.start),
                         inc := lv.inc));
                 AddSet(idxpool, asn.loc);
                 od;

             [asns, c] := Pull(c, [assign, @(0), [add, lv.v, @(2).cond(e->ObjId(e)<>Value)]],
                 e->skip(), e->e);
             for asn in asns do
                 Add(lvfifo, rec(v := asn.loc,
                         start := add(asn.exp.args[2], lv.start),
                         inc := lv.inc));
                 AddSet(idxpool, asn.loc);
                 od;

             [asns, c] := Pull(c, [assign, @(0), [mul, @(1), lv.v]],
                 e->skip(), e->e);
             for asn in asns do
                 Add(lvfifo, rec(v := asn.loc,
                     start := mul(asn.exp.args[1], lv.start),
                     inc := mul(asn.exp.args[1],  lv.inc)));
                 AddSet(idxpool, asn.loc);
                 od;

             [asns, c] := Pull(c, [assign, @(0), [mul, lv.v, @(2)]],
                 e->skip(), e->e);
             for asn in asns do
                 Add(lvfifo, rec(v := asn.loc,
                     start := mul(asn.exp.args[2], lv.start),
                     inc := mul(asn.exp.args[2],  lv.inc)));
                 AddSet(idxpool, asn.loc);
                 od;

         od;

         #Some of the lvs were only used by other lvs, let's kick
         #them out
         freevars := Set(c.free());
         lvlist := Filtered(lvlist, lv -> lv.v in freevars);

         #Now we partition the lvs into bins that have the same increment
         #We call it an lvfamily
         lvfamilies:=[];
         for lv in lvlist do
             fams := Filtered(lvfamilies, fam->fam.inc=lv.inc);
             if (Length(fams)>0) then
                 Add(fams[1].list, lv);
             else
                 Add(lvfamilies, rec(inc:=lv.inc, list:=[lv]));
             fi;
         od;

         #If the increment is an operation, hoist this operation
         for lv in lvfamilies do
             if (ObjId(lv.inc) in [add, mul]) then
                d := var.fresh_t("d", lv.inc.t);
                Add(c.cmds, assign(d, lv.inc));
                lv.inc := d;
             fi;
         od;

         #OK, so this is the tricky part.
         #we want to fuse lvfamilies together if it is possible.
         #it is only possible if members of the family are used sequentially in the
         #loop AND if we know the pattern

         var_usage := [];
         DoForAll(c.rChildren(), function(x) if IsAssign(x) then Append(var_usage, x.op_in()); fi; end);

         for j in [1..Length(lvfamilies)] do
           fam:=lvfamilies[j];
           members:=List(fam.list, x->x.v);
           member_usage := RemoveAdjacentDuplicates(Filtered(var_usage, x -> x in members));

           if (Length(Set(member_usage))=Length(member_usage) and Length(members)>1) then
               #ok so this is a good candidate, but can we match it?
               #our target is that they all are add( add(param, i*param), param)
               rootexpr := fam.list[1].start;

               if (ObjId(rootexpr)=add) then
                   d := rootexpr.args[2];
                   equal_second_member:=true;
                   for m in fam.list do
                   if ((ObjId(m.start)<>add) or (m.start.args[2]<>d)) then
                       equal_second_member:=false;
                   fi;
                   od;

                   if ((equal_second_member) and (rootexpr.args[1] in fam.list[2].start.args[1].pred)) then
                           step := Difference(fam.list[2].start.args[1].pred, Set([rootexpr.args[1]]));
                           if (Length(step)=1) then
                               pointer := 2;
                               cc := [];
                               for d in c.rChildren() do
                               if ((Length(fam.list)>=pointer) and (members[pointer] in d.op_in())) then
                                   pointer:=pointer+1;
                              Add(cc, assign(members[1], add(members[1], step[1])));
                               fi;
                               Add(cc, d);
                               od;
                               c := chain(cc);

                               d := sub(fam.inc, mul(V(Length(members)-1), step[1]));
                               s := var.fresh_t("inc", d.t);
                               Add(cc, assign(s, d));

                               c := chain(cc);
                               lvfamilies[j] := rec(inc := s,
                                   list:=[rec(v:=members[1],
                                           start:=rootexpr,
                                           inc:=lvfamilies[j].inc)]);
                               for d in [2..Length(members)] do
                               SubstBottomUp(c, members[d], e->members[1]);
                               od;
                           fi;
                   fi;
               fi;
           fi;
        od;

         #This is the hoister logic
         candidate_memory_acc := Set([]);
         confirmed_memory_acc := Set([]);
         hoist := [];
         body := [];
         epilogue := [];

         for asn in c.cmds do
             if IsAssign(asn) then
                 s := asn.op_in();
                 IntersectSet(s, idxpool);
                 if (Length(s)=0) then
                     if (ObjId(asn.exp)=deref) then
                         #cannot assume SSA on derefs
                         #so we mark it and leave it
                         #we'll fix that later
                         AddSet(candidate_memory_acc, asn.exp);
                     fi;
                     Add(hoist, asn);
                 else
                     if ((ObjId(asn.loc)=deref) and (asn.loc in candidate_memory_acc)) then
                         AddSet(confirmed_memory_acc, asn.loc);
                     fi;
                     Add(body, asn);
                     UniteSet(idxpool, asn.op_out());
                 fi;
             else
                 if (ObjId(asn)<> skip) then
                     Error("Non assigns cannot be handled");
                 fi;
             fi;
         od;

         #Plug back the induced vars
         newloop_idx := [];
         for lv in lvfamilies do
            if ((newloop_idx=[])and(Length(lv.list)=1)) then
                newloop_idx := lv.list[1].v;
                d := lv.list[1].start;
                s := var.fresh_t("init", d.t);
                Add(hoist, assign(s, d));
                lv.list[1].start := s;
                d := add(lv.list[1].start, mul(loop_range, lv.list[1].inc));
                s := var.fresh_t("ubound", d.t);
                Add(hoist, assign(s, d));
                newloop_range := s;
                newloop_start := lv.list[1].start;
                newloop_inc := lv.inc;
            else
                for d in lv.list do
                   Add(hoist, assign(d.v, d.start));
                   Add(body, assign(d.v, add(d.v, lv.inc)));
                od;
            fi;
         od;
         if (newloop_idx=[]) then
             newloop_idx := loop_idx;
             newloop_range := loop_range;
             newloop_start := V(0);
             newloop_inc := V(1);
         fi;
         Add(hoist, assign(newloop_idx, newloop_start));
         Add(body, assign(newloop_idx, add(newloop_idx, newloop_inc)));

         #This is the memory accumulator logic
         for d in confirmed_memory_acc do
             asns := Collect(hoist, [assign, @(1), [deref, d.loc]]);
             if Length(asns)<>1 then
                 Error("Multi assign in hoist????");
             else
                 SubstBottomUp(body, [deref, d.loc], e->asns[1].loc);
                 Add(epilogue, assign(d, asns[1].loc));
             fi;
         od;

         c := Compile.declareVars(chain(
                 chain(hoist),
                 doloop(newloop_idx, newloop_range, chain(body)),
                 chain(epilogue)));
    fi;
    return c;
end;


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

#P Parametrized permutation classes
#P --------------------------------
#P

## -----------------------------------------------------------------------------
#F L(<n>, <str>) - stride permutation
##
Class(L, PermClass, rec(
    def := (n,str) -> Checked(IsPosIntSym(n), IsPosIntSym(str),
        (not (IsInt(n) and IsInt(str)) or n mod str = 0), rec()),

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

    lambda := self >> let(
        n := self.params[1], str := self.params[2], i := Ind(n),
        Lambda(i, idiv(i, n/str) + str * imod(i, n/str))),

    transpose := self >> self.__bases__[1](self.params[1], self.params[1] / self.params[2]),
    isSymmetric := self >> (self.params[1] = self.params[2]^2) or (self.params[2] = 1) or (self.params[2] = self.params[1]),
    printlatex := self >> Print(" \\stride^{", self.params[1], "}_{",self.params[2],"} ")

));

## -----------------------------------------------------------------------------
#F Tr(<k>, <str>) - stride permutation L(k*str, str),
##
## This L equivalent is used to simplify rewriting, since k is not explicit in L,
## and if symbolics are used, to extract it we have to do simplification tricks,
## eg.  n/f1(n) = f2(n).
##
## Tr stands for Transpose.
##
Class(Tr, PermClass, rec(
    def := (k,str) -> Checked(IsPosIntSym(k), IsPosIntSym(str), rec()),

    lambda := self >> let(
        k := self.params[1], str := self.params[2], i := Ind(k*str),
        Lambda(i, idiv(i, k) + str * imod(i, k))),

    domain := self >> self.params[1] * self.params[2],
    range := self >> self.params[1] * self.params[2],

    dims := self >> [self.range(), self.domain()],

    transpose := self >> self.__bases__[1](self.params[2], self.params[1])
));

## -----------------------------------------------------------------------------
#F Z(N,k) : Cyclic shift by <k>, permutation {0, ..., n-1} -> {k,...,n-1, 0,...,k-1}
##
# if one redefines .lambda, then .transposed should be also handled in custom .lambda
Class(Z, PermClass, rec(
    abbrevs := [ n -> [n,1] ],
    def := (n,k) -> rec(),
    domain := self >> self.params[1],
    range  := self >> self.params[1],
    lambda := self >> let(n:=self.params[1], k:=self.params[2], i := Ind(n),
        Lambda(i, imod(i + k, n))),
    transpose := self >> self.__bases__[1](self.params[1], self.params[1] - self.params[2])
));

## -----------------------------------------------------------------------------
#F J(N) : NxN reverse identity, also known as reversal permutation (1,n)(2,n-1)...
##
Class(J, PermClass, rec(
    def := n -> rec(),
    lambda := self >> let(i := Ind(self.params[1]), Lambda(i,self.params[1]-i-1)),
    domain := self >> self.params[1],
    range  := self >> self.params[1],
    transpose := self >> self,
    isSymmetric := True
));

## -----------------------------------------------------------------------------
#F OddStride(<n>, <str>)
#F OS(<n>, <str>)        - odd stride permutation (i -> i * str mod n)
##
Class(OS, PermClass, rec(
    isCyclic := true,
    def := (n,str) -> Checked(IsPosIntSym(n), IsIntSym(str), n > 0,
        AnySyms(n,str) or (Gcd(EvalScalar(n),EvalScalar(str))=1), rec()),
    domain := self >> self.params[1],
    range  := self >> self.params[1],
    lambda := self >> let(i := Ind(self.params[1]),
        Lambda(i, imod(self.params[2] * i, self.params[1]))),
    transpose := self >> self.__bases__[1](self.params[1], self.params[2]^-1 mod self.params[1])
));
OddStride := OS;

Declare(gammaTensor);
## -----------------------------------------------------------------------------
#F CRT(<r>, <s>) - Chinese remainder theorem permutation function of size r*s
##
Class(CRT, PermClass, rec(
    isCyclic := true,

    abbrevs := [ (r, s) -> Checked(IsPosIntSym(r), IsPosIntSym(s), AnySyms(r,s) or Gcd(r,s)=1,
        let(alpha := 1/s mod r,
        beta  := 1/r mod s,
        [r, s, alpha, beta])) ],

    toGammaTensor := self >> gammaTensor(
        OddStride(self.params[1], self.params[3]),
        OddStride(self.params[2], self.params[4])),

    def := (r, s, alpha, beta) -> Checked(
        IsPosIntSym(r), IsPosIntSym(s), AnySyms(r, s) or Gcd(r,s)=1,
        rec()),

    domain := self >> self.params[1] * self.params[2],
    range  := self >> self.params[1] * self.params[2],

    lambda := self >> let(
        r := self.params[1], s := self.params[2], N := r*s,
        alpha := self.params[3], beta := self.params[4],
        aa := (1/s/alpha) mod r, bb := (1/r/beta) mod s,
        i := Ind(N),
        When(not self.transposed,
            Lambda(i, ((s*alpha*idiv(i, s)) + (r*beta*imod(i, s))) mod N),
            Lambda(i,   s*(imod(i*aa, r))   + (imod(i*bb, s))))),
));

Class(fCond, FuncClass, rec(
    def := (cond, f1, f2) -> Checked(
        f1.range() = f2.range(),
        f1.domain() = f2.domain(),
	rec()),

    domain := self >> self.params[1].domain(),
    range := self >> self.params[2].range(),

    lambda := self >> let(f1 := self.params[2], f2:= self.params[3], mycond := self.params[1],
        i := Ind(f1.domain()),
        Lambda(i, cond(mycond.lambda().at(i), f1.lambda().at(i), f2.lambda().at(i)))),

    transpose := self >> let(base := self.__bases__[1],
        base(self.params[1], self.params[2].transpose(), self.params[3].transpose()))
));

# ----------------------------------------------------------
# Bit and Digit reversals
# ----------------------------------------------------------

# Return x as a b-bit value, formed as a vector
_numToBits := (k, b) >> Reversed(List([1..b], i-> imod(idiv(k, 2^(i-1)), 2)));

# Return number corresponding to bit vector b
_bitsToNum := b -> let(
   l := Reversed(b),
   Sum([1..Length(l)], i -> 2^(i-1) * l[i]));

# Calculate base-r digit reversed value of k, where k is between 0 and n-1
_digitRev := function(k, n, r)
   local bitval, numbits, bitsperchunk, numchunks, res_bits, i, s;

   numbits      := Log2Int(n);
   bitsperchunk := Log2Int(r);
   numchunks    := numbits / bitsperchunk;

   bitval := _numToBits(k, numbits);

   res_bits := [];
   for i in Reversed([1..numchunks]) do
       s := bitval{([bitsperchunk*i-(bitsperchunk-1) .. bitsperchunk * i])};
       Append(res_bits, s);
   od;

   return _bitsToNum(res_bits);
end;

#F DR(<k, r>) - digit reversal permutation R^k_r
##
Class(DR, PermClass, rec(
    def := (k, r) -> Checked(IsInt(k), Is2Power(k), (r^(LogInt(k,r)) = k), rec()),

    domain := self >> self.params[1],
    range  := self >> self.params[1],
 # This is very slow because of how Gap evaluates functions
 #    lambda := self >> let(
 #         n := self.params[1],  base := self.params[2],  t := LogInt(n, base),
 #         rev := List([1..t], i-> fTensor(fId(base^(t-i)), L(base^i, base))),
 #         fCompose(Reversed(rev)).lambda()),

# This is faster but a little bit uglier.
    lambda := self >> let(
        i := Ind(self.params[1]),
        Lambda(i, _digitRev(i, self.params[1], self.params[2]))),

    transpose := self >> self
));

#F BR(<k>) - bit reversal permutation
##
Class(BR, PermClass, rec(
    def := k -> Checked(IsInt(k), Is2Power(k),  rec()),
    domain := self >> self.params[1],
    range  := self >> self.params[1],
    lambda := self >> DR(self.params[1], 2).lambda(),
    transpose := self >> self
));


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

Class(SMP_UnparseMixin, rec(
    includes := ["<include/threads.h>", "<include/smp2.h>"],

    # nothing to do here, threads must be started outside
    smp_fork := (self, o, i, is) >> Print(
        Blanks(i), "/* region must be executed concurrently by ", o.nthreads, " threads */ \n",
        Blanks(i), "{\n",
        self(o.cmd, i+is, is),
        Blanks(i), "}\n"
    ),

    smp_loop := meth(self, o, i, is)
        local v, lo, hi;
        v := o.var;
        lo := 0;
        hi := o.range-1;
        Print(Blanks(i), "{ /* begin parallel loop */\n");
        Print(Blanks(i+is), self.printf("int $1 = $2; \n", [o.tidvar, o.tidexp]));
        Print(Blanks(i+is), self.printf("for(int $1 = $2; $1 <= $3; $1 += $4) {\n", [v, o.tidvar + lo, hi, o.nthreads]));
        self(o.cmd,i+is+is,is);
        Print(Blanks(i+is), "}\n");
        Print(Blanks(i), "} /* end parallel loop */\n");
    end,

    threadId := (self, o, i, is) >> Print("tid")   # Error ??
));

Class(OpenMP_UnparseMixin, SMP_UnparseMixin, rec(
    includes := ["<omp.h>"],

    # start threads using 'omp parallel' pragma
    smp_fork := (self, o, i, is) >> Print(
        Blanks(i), "#pragma omp parallel num_threads(", o.nthreads, ")\n",
        Blanks(i), "{\n",
        self(o.cmd, i+is, is),
        Blanks(i), "}\n"
    ),

    threadId := (self,o,i,is) >> Print("omp_get_thread_num()"),
    barrier := (self,o,i, is) >> Print("#pragma omp barrier\n")
));

Class(OpenMP_UnparseMixin_ParFor, SMP_UnparseMixin, rec(
    includes := ["<omp.h>"],

    # start threads using 'omp parallel' pragma
    smp_fork := (self, o, i, is) >>
        Print(Blanks(i), "/* SMP fork */\n",
        Blanks(i), "{\n",
#   NOTE: why do I need to do that??
        self.opts.unparser(o.cmd,i+is,is),
#        self(o.cmd, i+is, is),
        Blanks(i), "}\n"
    ),

    threadId := (self,o,i,is) >> Print("omp_get_thread_num()"),
    barrier := (self,o,i, is) >> Print(Blanks(i), "/* SMP barrier */\n"),

    smp_loop := (self,o,i,is) >> let(v := o.var, lo := 0, hi := o.range,
            Print(Blanks(i), "#pragma omp parallel for schedule(static, ", Int((hi+1)/2), ")\n",
            Blanks(i), "for(int ", v, " = ", lo, "; ", v, " < ", hi, "; ", v, "++) {\n",
                Blanks(i + is), "int ", o.tidvar, " = ", v, ";\n",
#   NOTE: why do I need to do that??
                self.opts.unparser(o.cmd,i+is,is),
#                self(o.cmd,i+is,is),
                Blanks(i), "}\n")),
));


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

# -------------------------------------------------------------------------
# SSE bit-level instructions
# 128x1i: 128-way bit register
# -------------------------------------------------------------------------

#Class(exp_128x1,  rec(v := 128, computeType := self >> TVect(T_UInt(1), 128)));
#Class(cmd_128x1,  rec(v := 128, computeType := self >> TVect(T_UInt(1), 128)));

Class(sseb_shl, VecExp_128.binary());
Class(sseb_shr, VecExp_128.binary());
Class(sseb_xor, VecExp_128.binary());
Class(sseb_and, VecExp_128.binary());
Class(sseb_or , VecExp_128.binary());
Class(sseb_not, VecExp_128.binary());

Class(sseb_bcast, VecExp_128.binary());

Class(sseb_uload,      VecExp_128.unary());
Class(sseb_load_sd_64, VecExp_128.unary());
Class(sseb_load_sd_32, VecExp_128.unary());
Class(sseb_load_sd_16, VecExp_128.unary());
Class(sseb_load_lo_64, VecExp_128.binary());
Class(sseb_load_hi_64, VecExp_128.binary());

Class(sseb_ustore,      VecExpCommand.binary());
Class(sseb_store_sd_64, VecExpCommand.binary());
Class(sseb_store_sd_32, VecExpCommand.binary());
Class(sseb_store_sd_16, VecExpCommand.binary());

Class(sseb_store_lo_64, VecExpCommand.ternary());
Class(sseb_store_hi_64, VecExpCommand.ternary());

#Class(vunpacklo_16x32f,  vbinop_new, exp_16x32f, rec(semantic := (in1, in2, p) -> unpacklo(in1, in2, 16, 1)));
#Class(vunpackhi_16x32f,  vbinop_new, exp_16x32f, rec(semantic := (in1, in2, p) -> unpackhi(in1, in2, 16, 1)));
#Class(vunpacklo2_16x32f, vbinop_new, exp_16x32f, rec(semantic := (in1, in2, p) -> unpacklo(in1, in2, 16, 2)));
#Class(vunpackhi2_16x32f, vbinop_new, exp_16x32f, rec(semantic := (in1, in2, p) -> unpackhi(in1, in2, 16, 2)));
#Class(vunpacklo4_16x32f, vbinop_new, exp_16x32f, rec(semantic := (in1, in2, p) -> unpacklo(in1, in2, 16, 4)));
#Class(vunpackhi4_16x32f, vbinop_new, exp_16x32f, rec(semantic := (in1, in2, p) -> unpackhi(in1, in2, 16, 4)));

# 128-bit mask split up in 8 bit chunks
 _vstore_16x8_mask := nmany -> List([0..15], i -> "0x" :: HexStringInt(Cond(
    nmany >= (i+1)*8, 255,
    nmany < i*8, 0,
    (nmany mod 8) <> 0, let(b := nmany mod 8,
	(2^b - 1) * 2^(8-b)),  # upper b bits (of 8) must be 1s
    0
)));

Class(SSEB_128, SIMD_Intel, rec(
    countrec := rec(
        ops := [
            [add, sub],
	    [mul],
            [sseb_bcast], # shuffles
            [sseb_load_sd_64, sseb_load_sd_32, sseb_load_sd_16],
            [deref],
            Value      # Value without [] is a keyword in countOps !!
        ],
        printstrings := ["[adds]", "[mults]", "[vperms]", "[svldst]", "[vldst]", "[vval]"],
        type := "TVect",
        arithcost := (self, opcount) >> opcount[1]+opcount[2]
    ),

    includes := () -> ["<include/bitsse.h>"] :: _MM_MALLOC() :: _PMMINTRIN() :: _EMMINTRIN() :: _XMMINTRIN(),
    active   := true,
    info     := "SSE2 128 x 1-bit",
    ctype    := "__bit",
    instr    := [],
    v        := 128,
    t        := BitVector(128),
    bits     := 1,
    isFloat  := false,
    isFixedPoint := false,
    splopts := rec(),

    v_ones  := BitVector(128).one(),
    v_zeros := BitVector(128).zero(),
    val     := bits -> BitVector(128).value(bits),

    dupload  := (self, y, x) >> Checked(ObjId(x)=nth, # NOTE: is there a better way?
	let(base := x.loc,
	    ofs  := x.idx,
	    v    := self.v,
	    xvec := Cond(IsUnalignedPtrT(base.t), sseb_uload(base + idiv(ofs,v)*v),
		         vtref(self.t, base, idiv(ofs, v))),
	    assign(y, sseb_bcast(xvec, imod(ofs, v))))),

    svload := [ [ ], # load using subvecs of len 1
                [ ], # load using subvecs of len 2
                [ ], # load using subvecs of len 4
    ],

    svstore := [ [ ], # store using subvecs of len 1
                 [ ], # store using subvecs of len 2
                 [ ], # store using subvecs of len 4
    ],

    # keep the n lower scalars and zero the other ones
    mask := (self, c, n) >> Cond(
	n = 128, c,
        bin_and(c, self.val(Replicate(n, 1) :: Replicate(self.v - n, 0)))),

    loadCont := (self, n, y, yofs, x, xofs, xofs_align, opts) >> let(
#	_t := Checked(xofs=0, 0), # NOTE: remove this limitation, this is because xofs is in bits
	a  := _unwrap(xofs_align),
	nn := _unwrap(n),
	yy := vtref(self.t, y, yofs),
	m  := x -> self.mask(x, nn),
	# NOTE: BELOW CODE IS NOT VERIFIED, in particular shift directions might be invalid
	# NOTE: sseb_shl/shr need to be implemented in the backend,
        # NOTE : SSE requires shift values to be immediates (or does it not?)
	Cond(a = 0 and not IsUnalignedPtrT(x.t),
	         assign(yy, m(vtref(self.t, x, div(xofs, self.v)))),

	     # known alignment, sv is small, so that we only need 1 aligned load + 1 shift + mask
	     IsInt(a) and not IsUnalignedPtrT(x.t) and (nn <= self.v - a),
		 let(v1 := vtref(self.t, x, div(xofs - a, self.v)),
		     assign(yy, m(sseb_shl(v1, a)))),

	     # known alignment, sv covers 2 vectors, 2 aligned loads + 2 shifts + mask
	     # NB: no masking is needed because shifts will do the job
	     IsInt(a) and not IsUnalignedPtrT(x.t),
		 let(v1 := vtref(self.t, x, div(xofs - a, self.v)),
		     v2 := vtref(self.t, x, div(xofs - a, self.v) + 1),
		     assign(yy, m(sseb_or(sseb_shl(v1, a), sseb_shr(v2, self.v - a))))),

             # else, unknown alignment, use unaligned load
	     assign(yy, m(sseb_uload(x + xofs)))
        )
    ),

    # NOTE: The store mask has 8-bit granularity, how do we express this??
    storeCont := (self, n, y, yofs, yofs_align, x, xofs, opts) >> let(
#	_t := Checked(yofs=0, 0), # NOTE: remove this limitation, this is because yofs is in bits
	a  := _unwrap(yofs_align),
	nn := _unwrap(n),
	xx := vtref(self.t, x, xofs),
	Cond(nn = self.v and a = 0 and not IsUnalignedPtrT(y.t),
	        assign(vtref(self.t, y, div(yofs, self.v)), xx),
	     nn = self.v,
	        sseb_ustore(y + yofs, xx),
	     # else
	     #NOTE: yofs is in bits:
	        vstoremsk_16x8i(nth(y, yofs).toPtr(TReal), xx, _vstore_16x8_mask(nn)))),

));

_bitsBytes := (bits) -> let(bytelsts := List([0..Length(bits)/8-1], i -> bits{ [1..8] + 8*i }),
	                 List(bytelsts, l -> Sum([ 1 .. 8 ], i -> l[i] * 2 ^ (i - 1))));

Class(BitSSEUnparser, SSEUnparser, rec(
    infixbreak := 20,

    BitVector := (self, t, vars, i, is) >> Print("__m128i ", self.infix(vars, ", ", i + is)),

    Value := (self, o, i, is) >> Cond(ObjId(o.t) = BitVector,
	let(bytenums := _bitsBytes(o.v),
            Cond(self.cx.isInside(Value) and Length(self.cx.Value) >= 2, # nested in an array
		 Print(            "{", self.infix(bytenums, ", "), "}"),
		 Print("_mm_set_epi8(", self.infix(Reversed(bytenums), ", "), ")"))),

	Inherited(o, i, is)),

    mul := (self, o, i, is) >> Cond(ObjId(o.t) = BitVector and Length(o.args)=2,
	self.prefix("_mm_and_si128", o.args),
	Inherited(o, i, is)),

    add := (self, o, i, is) >> Cond(ObjId(o.t) = BitVector and Length(o.args)=2,
	self.prefix("_mm_xor_si128", o.args),
	Inherited(o, i, is)),

    sseb_uload := (self, o, i, is) >> self(vloadu_16x8i(o.args[1]), i, is),

    sseb_ustore := (self, o, i, is) >> self(vstoreu_16x8i(o.args[1], o.args[2]), i, is)
));


text:
Class(GathIntoVec, BaseOperation, rec(
  __call__ := (self, sr, var, range) >> SPL(WithBases(self, rec(
    sr := Checked(IsSemiring(sr), sr),
	var := var,
	range := range,
    _children := [range]
  ))),
  dims := self >> [self.var, V(1)],
  from_rChildren := (self, rch) >> CopyFields(self, rec(_children := rch)),
  doNotMarkBB := true,
  print := (self, i, is) >> Print(
        self.name, "(", self.sr, ", ", self.var, ", ", self.range, ")"
  )
));

DefaultSumsGen.GathIntoVec := (self, o, opts) >> o;

DefaultCodegen.GathIntoVec := meth(self, o, y, x, opts)
  local i, x_first, x_last;
  if (IsBound(opts.isCSR) and opts.isCSR) then
    i := var.fresh_t("i", x.t);
    x_first := var.fresh_t("x_first", x.t);
    x_last := var.fresh_t("x_last", x.t);
	return	decl([i, x_first, x_last], chain (
		assign(x_first, add(x, o.var)),
		assign(x_last, add(x, add(o.var, V(1)))),

		loopf(i, deref(x_first), deref(x_last)),
		  chain(
			assign(deref(add(y,sub(i,x_first))), sub(deref(x_last), deref(x_first)))
		  )
		)
	);
  fi;
end;