text:
#
# Read("~/Workspace/groupsSB/epi/group.gi");
#
#
# needs:
# type:="A";
# rank:=2;
# nr_pos_roots:=3;
#

ZZ:=Integers;
avarnames:=List([1..100],i->Concatenation("a_{",String(i),"}"));
bvarnames:=List([1..100],i->Concatenation("b_{",String(i),"}"));
cvarnames:=List([1..100],i->Concatenation("c_{",String(i),"}"));
xvarnames:=List([1..100],i->Concatenation("x_{",String(i),"}"));
varnames:=Concatenation(avarnames,bvarnames,cvarnames,xvarnames);
APR:=PolynomialRing(ZZ,varnames);
vars:=IndeterminatesOfPolynomialRing(APR);
xvars:=vars{[301..400]};

sla:=SimpleLieAlgebraTypeA_G(type,rank,APR);

cb:=CanonicalBasis(sla);

e:=cb[1];

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

ade:=function(e)
	local result,v;
	result:=[];
	for v in cb do
		Append(result,[Coefficients(cb,e*v)]);
	od;
	result:=TransposedMat(result);
	return result;
end;

root_group:=function(index,t)
	local ee,tmp,result,i;
	ee:=ade(cb[index]);
	tmp:=ee;
	result:=One(APR)*tmp^0;
	i:=1;
	while Length(Set(Concatenation(tmp)))<>1 do
		result:=result+t^i*tmp/Factorial(i);
		i:=i+1;
		tmp:=tmp*ee;
	od;
	return result;
end;
#u1a1:=root_group(1,vars[1]);

pos_root_groups:=function(start_a_index)
	return List([1..nr_pos_roots],i->root_group(i,vars[start_a_index+i]));
end;

generic_U:=function(start_a_index)
	local Uas;
	Uas:=pos_root_groups(start_a_index);
	return Product(Uas);
end;

Ua:=generic_U(0);
Ub:=generic_U(10);
Uc:=generic_U(20);

#
#
#

evaluate_U:=function(u,vals)
	local i,j,result,v;
	result := [];
	for i in [1..Length(u)] do
		Append(result,[[1..Length(u)]]);
		for j in [1..Length(u)] do
			result[i][j]:=u[i][j];
		od;
	od;
	Print(result);
	for i in [1..Length(u)] do
		for j in [1..Length(u)] do
			for v in vals do
				result[i][j]:=One(APR)*Value(One(APR)*result[i][j],v[1],v[2]);
			od;
		od;
	od;
	return result;
end;

evaluate_U:=function(u,vals)
	local i,j,result,v;
	result := [];
	for i in [1..Length(u)] do
		Append(result,[[1..Length(u)]]);
		for j in [1..Length(u)] do
			result[i][j]:=u[i][j];
		od;
	od;
	Print(result);
	for i in [1..Length(u)] do
		for j in [1..Length(u)] do
			for v in vals do
				result[i][j]:=One(APR)*Value(One(APR)*result[i][j],v[1],v[2]);
			od;
		od;
	od;
	return result;
end;

evaluate_rels:=function(rels,vals)
	local i,result,v;
	result :=List([1..Length(rels)],i->rels[i]);
	for i in [1..Length(rels)] do
		for v in vals do
			#Print(Length(rels),": ",rels[i],"\n");
			result[i]:=One(APR)*Value(result[i],v[1],v[2]);
			#nn[i]:=One(APR)*Value(nn[i],v[1],v[2]);
		od;
	od;
	return result;
end;


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

benchNEON := function()
    return rec(
        2x32f := rec(
            wht := rec(
                small := _defaultSizes(s->doSimdWht(s, NEON_HALF, rec(
		            verify := true, oddSizes := true, svct := true, stdTTensor := true, tsplPFA := false)), [4]),
                medium := _defaultSizes(s->doSimdWht(s, NEON_HALF, rec(
			    oddSizes := false, svct := true, stdTTensor := true, tsplPFA := false)), List([2..10], i->2^i))
                ),
            1d := rec(
                dft_sc := rec(
                    small := _defaultSizes(s->doSimdDft(s, NEON_HALF, rec(
				verify:=true, tsplBluestein:=false, interleavedComplex := false, PRDFT:=true, URDFT:= true,
				cplxVect := true, stdTTensor := false, globalUnrolling:=10000)),
                        [ 2..32 ]),
                    medium := _defaultSizes(s->doSimdDft(s, NEON_HALF, rec(tsplRader:=false, tsplBluestein:=false,
				tsplPFA:=false, oddSizes:=false, interleavedComplex := false)),
                        _svctSizes(1024, 16, 2)),

                    ),
                dft_ic := rec(
                    small := _defaultSizes(s->doSimdDft(s, NEON_HALF, rec(
				verify:=true, tsplBluestein:=false, interleavedComplex := true, PRDFT:=true, URDFT:= true,
				cplxVect := true, stdTTensor := false, globalUnrolling:=10000)),
                        [2..32]),
                    medium := _defaultSizes(s->doSimdDft(s, NEON_HALF, rec(
				tsplRader:=false, tsplBluestein:=false, tsplPFA:=false, oddSizes:=false, interleavedComplex := true)),
                        _svctSizes(1024, 16, 2)),
                    medium_cx := _defaultSizes(s->doSimdDft(s, NEON_HALF, rec(
				tsplRader:=false, tsplBluestein:=false, tsplPFA:=false, oddSizes:=false, interleavedComplex := true,
				cplxVect := true, realVect := false)),
                        _svctSizes(1024, 16, 2)),
                    ),

                trdft := _defaultSizes(s->doSimdSymDFT(TRDFT, s, NEON_HALF, rec( verify:=true,
                            interleavedComplex := true, PRDFT:=true, URDFT:= true, tsplBluestein := false, cplxVect := true,
                            realVect := true, propagateNth := true, useDeref := true,
                            globalUnrolling:=10000)), 4*[1..32]),

                dht := _defaultSizes(s->doSimdSymDFT(TDHT, s, NEON_HALF, rec(verify :=true)), 8*[1..12]),
                dct2 := _defaultSizes(s->doSimdSymDFT(TDCT2, s, NEON_HALF, rec(verify :=true)), 8*[1..12]),
                dct3 := _defaultSizes(s->doSimdSymDFT(TDCT3, s, NEON_HALF, rec(verify :=true)), 8*[1..12]),
                dct4 := _defaultSizes(s->doSimdSymDFT(TDCT4, s, NEON_HALF, rec(verify :=true)), 8*[1..12]),
                dst2 := _defaultSizes(s->doSimdSymDFT(TDST2, s, NEON_HALF, rec(verify :=true)), 8*[1..12]),
                dst3 := _defaultSizes(s->doSimdSymDFT(TDST3, s, NEON_HALF, rec(verify :=true)), 8*[1..12]),
                dst4 := _defaultSizes(s->doSimdSymDFT(TDST4, s, NEON_HALF, rec(verify :=true)), 8*[1..12]),
                mdct := _defaultSizes(s->doSimdSymDFT(TMDCT, s, NEON_HALF, rec(verify :=true)), 8*[1..12]),
                imdct := _defaultSizes(s->doSimdSymDFT(TIMDCT, s, NEON_HALF, rec(verify :=true)), 8*[1..12])
                ),

            2d := rec(
                dft_ic := rec(
                    medium := _defaultSizes(s -> doSimdMddft(s, NEON_HALF, rec(
				interleavedComplex := true,
                                oddSizes := false, svct := true, splitL := false, pushTag := true,
				flipIxA := false, stdTTensor := true, tsplPFA := false)),
                        4 * List([1..16], i->[i,i])),

                    small := _defaultSizes(s->doSimdMddft(s, NEON_HALF, rec(
				verify:=true, interleavedComplex := true, globalUnrolling:=10000,
                                tsplPFA := false, pushTag:= false, oddSizes := true, svct := true, splitL := false)),
                        List([2..16], i->[i,i]))
                    ),
                dft_sc := rec(
                    medium := _defaultSizes(s -> doSimdMddft(s, NEON_HALF, rec(
				interleavedComplex := false,
                                oddSizes := false, svct := true, splitL := false, pushTag := true, flipIxA := false,
				stdTTensor := true, tsplPFA := false)),
                        4*List([1..16], i->[i,i])),
                    small := _defaultSizes(s->doSimdMddft(s, NEON_HALF, rec(verify:=true, interleavedComplex := false,
				globalUnrolling:=10000, tsplPFA := false, pushTag:= false, oddSizes := true,
				svct := true, splitL := false)),
                        List([2..16], i->[i,i]))
                    ),
                dct2 := _defaultSizes(s->doSimdSymMDDFT(DCT2, s, NEON_HALF, rec(verify := true)), [4, 8, 12, 16]),
                dct3 := _defaultSizes(s->doSimdSymMDDFT(DCT3, s, NEON_HALF, rec(verify := true)), [4, 8, 12, 16]),
                dct4 := _defaultSizes(s->doSimdSymMDDFT(DCT4, s, NEON_HALF, rec(verify := true)), [4, 8, 12, 16]),
                dst2 := _defaultSizes(s->doSimdSymMDDFT(DST2, s, NEON_HALF, rec(verify := true)), [4, 8, 12, 16]),
                dst3 := _defaultSizes(s->doSimdSymMDDFT(DST3, s, NEON_HALF, rec(verify := true)), [4, 8, 12, 16]),
                dst4 := _defaultSizes(s->doSimdSymMDDFT(DST4, s, NEON_HALF, rec(verify := true)), [4, 8, 12, 16])
            )
            ),

        4x32f := rec(
            wht := rec(
                small := _defaultSizes(s -> doSimdWht(s, NEON, rec(
			    verify := true, oddSizes := true, svct := true, stdTTensor := true, tsplPFA := false)),
		    List([1..3], i->2^i)),
                medium := _defaultSizes(s -> doSimdWht(s, NEON, rec(
			    oddSizes := false, svct := true, stdTTensor := true, tsplPFA := false
			    )), List([4..10], i->2^i))
                ),
            1d := rec(
                dft_sc := rec(
                    small := _defaultSizes(s->doSimdDft(s, NEON, rec(
				verify:=true, interleavedComplex := false, stdTTensor := false, globalUnrolling:=10000)),
                        [ 2..64 ]),
                    medium := _defaultSizes(s->doSimdDft(s, NEON, rec(
				tsplRader:=false, tsplBluestein:=false, tsplPFA:=false, oddSizes:=false, interleavedComplex := false)),
                        _svctSizes(1024, 16, 4)),

                    ),
                dft_ic := rec(
                    small := _defaultSizes(s->doSimdDft(s, NEON, rec(
				verify:=true, interleavedComplex := true, PRDFT:=true, URDFT:= true,
				cplxVect := true, stdTTensor := false, globalUnrolling:=10000)),
                        [ 2..64 ]),
                    medium := _defaultSizes(s->doSimdDft(s, NEON, rec(
				tsplRader:=false, tsplBluestein:=false, tsplPFA:=false,
				oddSizes:=false, interleavedComplex := true)),
                        _svctSizes(1024, 16, 4)),
                    medium_cx := _defaultSizes(s->doSimdDft(s, NEON, rec(
				tsplRader:=false, tsplBluestein:=false, tsplPFA:=false,
				oddSizes:=false, interleavedComplex := true,
                                cplxVect := true, realVect := false, PRDFT := false, URDFT := true)),
                        _svctSizes(1024, 16, 4))
                    ),

                rdft := _defaultSizes(s->doSimdSymDFT(TRDFT, s, NEON, rec(verify :=true)), 32*[1..12]),

                trdft := _defaultSizes(s->doSimdSymDFT(TRDFT, s, NEON, rec(
			    verify := true, interleavedComplex := true, PRDFT:=true, URDFT:= true,
			    tsplBluestein := false, cplxVect := true,   stdTTensor := false,
			    realVect := true, propagateNth := true, useDeref := true, globalUnrolling:=10000)),
		    8*[1..32]),

                dht := _defaultSizes(s->doSimdSymDFT(TDHT, s, NEON, rec(verify :=true)), 32*[1..12]),
                dct2 := _defaultSizes(s->doSimdSymDFT(TDCT2, s, NEON, rec(verify :=true)), 32*[1..12]),
                dct3 := _defaultSizes(s->doSimdSymDFT(TDCT3, s, NEON, rec(verify :=true)), 32*[1..12]),
                dct4 := _defaultSizes(s->doSimdSymDFT(TDCT4, s, NEON, rec(verify :=true)), 32*[1..12]),
                dst2 := _defaultSizes(s->doSimdSymDFT(TDST2, s, NEON, rec(verify :=true)), 32*[1..12]),
                dst3 := _defaultSizes(s->doSimdSymDFT(TDST3, s, NEON, rec(verify :=true)), 32*[1..12]),
                dst4 := _defaultSizes(s->doSimdSymDFT(TDST4, s, NEON, rec(verify :=true)), 32*[1..12]),
                mdct := _defaultSizes(s->doSimdSymDFT(TMDCT, s, NEON, rec(verify :=true)), 32*[1..12]),
                imdct := _defaultSizes(s->doSimdSymDFT(TIMDCT, s, NEON, rec(verify :=true)), 32*[1..12])
                ),
            2d := rec(
                dft_ic := rec(
                    medium := _defaultSizes(s -> doSimdMddft(s, NEON, rec(
				interleavedComplex := true,
                                oddSizes := false, svct := true, splitL := false, pushTag := true,
				flipIxA := false, stdTTensor := true, tsplPFA := false)),
                        16*List([1..8], i->[i,i])),
                    small := _defaultSizes(s->doSimdMddft(s, NEON, rec(
				verify:=true, interleavedComplex := true, globalUnrolling:=10000,
                                tsplPFA := false, pushTag:= false, oddSizes := true, svct := true, splitL := false)),
                        List([2..16], i->[i,i]))
                    ),
                dft_sc := rec(
                    medium := _defaultSizes(s -> doSimdMddft(s, NEON, rec(
				interleavedComplex := false,
                                oddSizes := false, svct := true, splitL := false, pushTag := true,
				flipIxA := false, stdTTensor := true, tsplPFA := false)),
                        16 * List([1..8], i->[i,i])),
                    small := _defaultSizes(s->doSimdMddft(s, NEON, rec(verify:=true, interleavedComplex := false,
				globalUnrolling:=10000, tsplPFA := false, pushTag:= false, oddSizes := true,
				svct := true, splitL := false)),
                        List([2..16], i->[i,i]))
                    ),
                dct2 := _defaultSizes(s -> doSimdSymMDDFT(DCT2, s, NEON, rec(verify := true)), [4, 8, 12, 16]),
                dct3 := _defaultSizes(s -> doSimdSymMDDFT(DCT3, s, NEON, rec(verify := true)), [4, 8, 12, 16]),
                dct4 := _defaultSizes(s -> doSimdSymMDDFT(DCT4, s, NEON, rec(verify := true)), [4, 8, 12, 16]),
                dst2 := _defaultSizes(s -> doSimdSymMDDFT(DST2, s, NEON, rec(verify := true)), [4, 8, 12, 16]),
                dst3 := _defaultSizes(s -> doSimdSymMDDFT(DST3, s, NEON, rec(verify := true)), [4, 8, 12, 16]),
                dst4 := _defaultSizes(s -> doSimdSymMDDFT(DST4, s, NEON, rec(verify := true)), [4, 8, 12, 16]),
                conv := _defaultSizes(s -> let(
                        opts := SIMDGlobals.getOpts(NEON, rec(svct := true, oddSizes := false,
				splitComplexTPrm := true, TRCDiag_VRCLR := true, globalUnrolling := 150,
				measureFinal := false)),
                        t := Flat(List(s, n -> let(
				    tt := TRConv2D(ImageVar([n, n])).withTags(opts.tags), [tt, tt.forwardTransform()]))),
                        spiral.libgen.DPBench.build(t, opts)), 16*[1..32])
            )
        )
    );
end;

ext: gi


text:
# Built-in
IdentityMat(3);

# One can also specify the base ring
IdentityMat(3, Integers mod 10);

ext: gap

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
#
#
#

Read("~/Workspace/Chevalley.gap/init.gi");

Read(Filename(home_dir,"lib/poly.gd"));
Read(Filename(home_dir,"lib/poly.gi"));


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

#---------------------------------------------------------------------------------------
#   SPU vector instructions
#---------------------------------------------------------------------------------------

Class(vbinop_8x16i_spu, vbinop_spu, rec( v := self >> 8));
Class(vbinop_4x32f_spu, vbinop_spu, rec( v := self >> 4));
Class(vbinop_2x64f_spu, vbinop_spu, rec( v := self >> 2));

Class(exp_8x16i, rec(v := 8, computeType := self >> TVect(TInt,    8)));
Class(exp_4x32f, rec(v := 4, computeType := self >> TVect(TDouble, 4)));
Class(exp_2x64f, rec(v := 2, computeType := self >> TVect(TDouble, 2)));

Class(cmd_8x16i, rec(v := 8, computeType := self >> TVect(TInt,    8)));
Class(cmd_4x32f, rec(v := 4, computeType := self >> TVect(TDouble, 4)));
Class(cmd_2x64f, rec(v := 2, computeType := self >> TVect(TDouble, 2)));

# Load -----------------------------
Class(vloadu8_spu8x16i, vloadop_new, exp_8x16i, rec(numargs := 1));
Class(vloadu4_spu4x32f, vloadop_new, exp_4x32f, rec(numargs := 1));

# Store ----------------------------

# Zero -----------------------------
Class(vzero_8x16i, vop_new, exp_8x16i, rec(numargs := 0));
Class(vzero_4x32f, vop_new, exp_4x32f, rec(numargs := 0));
Class(vzero_2x64f, vop_new, exp_2x64f, rec(numargs := 0));

# Subvec ---------------------------
Class(promote_spu8x16i, vbinop_new, exp_8x16i);
Class(promote_spu4x32f, vbinop_new, exp_4x32f);
Class(promote_spu2x64f, vbinop_new, exp_2x64f);

Class(extract_spu8x16i, vbinop_new, exp_8x16i, rec(computeType := self >> TInt));
Class(extract_spu4x32f, vbinop_new, exp_4x32f, rec(computeType := self >> TReal));
Class(extract_spu2x64f, vbinop_new, exp_2x64f, rec(computeType := self >> TReal));

Class(insert_spu8x16i, vop_new, exp_8x16i, rec(numargs := 3));
Class(insert_spu4x32f, vop_new, exp_4x32f, rec(numargs := 3));
Class(insert_spu2x64f, vop_new, exp_2x64f, rec(numargs := 3));

Class(vsplat_8x16i, vloadop_new, exp_8x16i, rec(numargs := 1));
Class(vsplat_4x32f, vloadop_new, exp_4x32f, rec(numargs := 1));
Class(vsplat_2x64f, vloadop_new, exp_2x64f, rec(numargs := 1));

# Binary ---------------------------
# VA: This breaks bin_or, looks like something unfinished so commented out.
# Class(bin_or, vbinop_new, exp_8x16i);
# Class(bin_or, vbinop_new, exp_4x32f);
# Class(bin_or, vbinop_new, exp_2x64f);

# Rotates ---------------------------
Class(slqwbyte_spu4x32f, vop_new, exp_4x32f, rec(numargs := 1));

Class(rlmaskqwbyte_spu4x32f, vop_new, exp_4x32f, rec(numargs := 1));

# Binary shuffle -------------------
#NOTE: What are sparams, params, semantic, and permparams?
#NOTE: Can we combine all these perms into the same thing somehow?
Class(vperm_8x16i_spu, vbinop_8x16i_spu, rec(
    semantic := (in1, in2, p) -> vpermop(in1, in2, p, 8),
    params := self >> sparams(8, 16),
    from_rChildren := (self, rch) >> ApplyFunc(ObjId(self), [rch[1], rch[2], rch[3].p]),
    permparams := aperm
));

Class(vperm_4x32f_spu, vbinop_4x32f_spu, rec(
    semantic := (in1, in2, p) -> vpermop(in1, in2, p, 4),
    params := self >> sparams(4, 8),

    #HACK: small hack: it'd be nice to not define from_rChildren explicitly
    #here. We have to do it though, because the last param is a perm
    #(vparam_spu) type, but the object must be created with a List, and not a
    #vparam_spu.

    from_rChildren := (self, rch) >> ApplyFunc(ObjId(self), [rch[1], rch[2], rch[3].p]),
    permparams := aperm
));

Class(vperm_2x64f_spu, vbinop_2x64f_spu, rec(
    semantic := (in1, in2, p) -> vpermop(in1, in2, p, 2),
    params := self >> sparams(2, 4),
    from_rChildren := (self, rch) >> ApplyFunc(ObjId(self), [rch[1], rch[2], rch[3].p]),
    permparams := aperm
));

# Unary shuffle --------------------
Class(vuperm_8x16i_spu, vunbinop_spu,
    rec(binop := vperm_8x16i_spu,
    from_rChildren := (self, rch) >> ApplyFunc(ObjId(self), [rch[1], rch[2].p])
));

Class(vuperm_4x32f_spu, vunbinop_spu,
    rec(binop := vperm_8x16i_spu,
    from_rChildren := (self, rch) >> ApplyFunc(ObjId(self), [rch[1], rch[2].p])
));

Class(vuperm_2x64f_spu, vunbinop_spu,
    rec(binop := vperm_8x16i_spu,
    from_rChildren := (self, rch) >> ApplyFunc(ObjId(self), [rch[1], rch[2].p])
));


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

SIMD_ISA_DB.hashAdd(TL(4, 2, 1, 2).withTags([ AVecReg(spu_4x32f) ]), [ rec(
      ruletree := SIMD_ISA_Bases1( TL(4, 2, 1, 2).withTags([ AVecReg(ppu_4x32f) ]) ),
      origtree := SIMD_ISA_Bases1( TL(4, 2, 1, 2).withTags([ AVecReg(ppu_4x32f) ]) ),
      measured := 0,
      globalUnrolling := 8 ) ]);
SIMD_ISA_DB.hashAdd(TL(4, 2, 2, 1).withTags([ AVecReg(spu_4x32f) ]), [ rec(
      ruletree := SIMD_ISA_Bases1( TL(4, 2, 2, 1).withTags([ AVecReg(ppu_4x32f) ]) ),
      origtree := SIMD_ISA_Bases1( TL(4, 2, 2, 1).withTags([ AVecReg(ppu_4x32f) ]) ),
      measured := 0,
      globalUnrolling := 8 ) ]);
SIMD_ISA_DB.hashAdd(TL(4, 2, 1, 4).withTags([ AVecReg(spu_4x32f) ]), [ rec(
      ruletree := IxLxI_vtensor( TL(4, 2, 1, 4).withTags([ AVecReg(ppu_4x32f) ]) ),
      origtree := IxLxI_vtensor( TL(4, 2, 1, 4).withTags([ AVecReg(ppu_4x32f) ]) ),
      measured := 0,
      globalUnrolling := 8 ) ]);
SIMD_ISA_DB.hashAdd(TL(4, 2, 4, 1).withTags([ AVecReg(spu_4x32f) ]), [ rec(
      ruletree := SIMD_ISA_Bases1( TL(4, 2, 4, 1).withTags([ AVecReg(ppu_4x32f) ]) ),
      origtree := SIMD_ISA_Bases1( TL(4, 2, 4, 1).withTags([ AVecReg(ppu_4x32f) ]) ),
      measured := 0,
      globalUnrolling := 8 ) ]);
SIMD_ISA_DB.hashAdd(TL(16, 8, 1, 1).withTags([ AVecReg(spu_4x32f) ]), [ rec(
      ruletree := IxLxI_kmn_km( TL(16, 8, 1, 1).withTags([ AVecReg(ppu_4x32f) ]),
          SIMD_ISA_Bases1( TL(4, 2, 4, 1).withTags([ AVecReg(ppu_4x32f) ]) ),
          SIMD_ISA_Bases2( TL(8, 4, 1, 2).withTags([ AVecReg(ppu_4x32f) ]) ) ),
      origtree := IxLxI_kmn_km( TL(16, 8, 1, 1).withTags([ AVecReg(ppu_4x32f) ]),
          SIMD_ISA_Bases1( TL(4, 2, 4, 1).withTags([ AVecReg(ppu_4x32f) ]) ),
          SIMD_ISA_Bases2( TL(8, 4, 1, 2).withTags([ AVecReg(ppu_4x32f) ]) ) ),
      measured := 0,
      globalUnrolling := 10000 ) ]);
SIMD_ISA_DB.hashAdd(TL(8, 4, 1, 1).withTags([ AVecReg(spu_4x32f) ]), [ rec(
      ruletree := SIMD_ISA_Bases1( TL(8, 4, 1, 1).withTags([ AVecReg(ppu_4x32f) ]) ),
      origtree := SIMD_ISA_Bases1( TL(8, 4, 1, 1).withTags([ AVecReg(ppu_4x32f) ]) ),
      measured := 0,
      globalUnrolling := 8 ) ]);
SIMD_ISA_DB.hashAdd(TL(8, 4, 1, 2).withTags([ AVecReg(spu_4x32f) ]), [ rec(
      ruletree := SIMD_ISA_Bases2( TL(8, 4, 1, 2).withTags([ AVecReg(ppu_4x32f) ]) ),
      origtree := SIMD_ISA_Bases2( TL(8, 4, 1, 2).withTags([ AVecReg(ppu_4x32f) ]) ),
      measured := 0,
      globalUnrolling := 8 ) ]);
SIMD_ISA_DB.hashAdd(TL(8, 4, 2, 1).withTags([ AVecReg(spu_4x32f) ]), [ rec(
      ruletree := SIMD_ISA_Bases1( TL(8, 4, 2, 1).withTags([ AVecReg(ppu_4x32f) ]) ),
      origtree := SIMD_ISA_Bases1( TL(8, 4, 2, 1).withTags([ AVecReg(ppu_4x32f) ]) ),
      measured := 0,
      globalUnrolling := 8 ) ]);
SIMD_ISA_DB.hashAdd(TL(8, 2, 1, 2).withTags([ AVecReg(spu_4x32f) ]), [ rec(
      ruletree := SIMD_ISA_Bases2( TL(8, 2, 1, 2).withTags([ AVecReg(ppu_4x32f) ]) ),
      origtree := SIMD_ISA_Bases2( TL(8, 2, 1, 2).withTags([ AVecReg(ppu_4x32f) ]) ),
      measured := 0,
      globalUnrolling := 8 ) ]);
SIMD_ISA_DB.hashAdd(TL(8, 2, 2, 1).withTags([ AVecReg(spu_4x32f) ]), [ rec(
      ruletree := SIMD_ISA_Bases1( TL(8, 2, 2, 1).withTags([ AVecReg(ppu_4x32f) ]) ),
      origtree := SIMD_ISA_Bases1( TL(8, 2, 2, 1).withTags([ AVecReg(ppu_4x32f) ]) ),
      measured := 0,
      globalUnrolling := 8 ) ]);
SIMD_ISA_DB.hashAdd(TL(16, 4, 1, 1).withTags([ AVecReg(spu_4x32f) ]), [ rec(
      ruletree := IxLxI_kmn_n( TL(16, 4, 1, 1).withTags([ AVecReg(ppu_4x32f) ]),
          SIMD_ISA_Bases2( TL(8, 4, 1, 2).withTags([ AVecReg(ppu_4x32f) ]) ),
          SIMD_ISA_Bases1( TL(8, 4, 2, 1).withTags([ AVecReg(ppu_4x32f) ]) ) ),
      origtree := IxLxI_kmn_n( TL(16, 4, 1, 1).withTags([ AVecReg(ppu_4x32f) ]),
          SIMD_ISA_Bases2( TL(8, 4, 1, 2).withTags([ AVecReg(ppu_4x32f) ]) ),
          SIMD_ISA_Bases1( TL(8, 4, 2, 1).withTags([ AVecReg(ppu_4x32f) ]) ) ),
      measured := 0,
      globalUnrolling := 8 ) ]);
SIMD_ISA_DB.hashAdd(TL(4, 2, 2, 2).withTags([ AVecReg(spu_4x32f) ]), [ rec(
      ruletree := SIMD_ISA_Bases1( TL(4, 2, 2, 2).withTags([ AVecReg(ppu_4x32f) ]) ),
      origtree := SIMD_ISA_Bases1( TL(4, 2, 2, 2).withTags([ AVecReg(ppu_4x32f) ]) ),
      measured := 0,
      globalUnrolling := 8 ) ]);
SIMD_ISA_DB.hashAdd(TL(16, 2, 1, 1).withTags([ AVecReg(spu_4x32f) ]), [ rec(
      ruletree := IxLxI_kmn_n( TL(16, 2, 1, 1).withTags([ AVecReg(ppu_4x32f) ]),
          SIMD_ISA_Bases2( TL(8, 2, 1, 2).withTags([ AVecReg(ppu_4x32f) ]) ),
          SIMD_ISA_Bases1( TL(4, 2, 4, 1).withTags([ AVecReg(ppu_4x32f) ]) ) ),
      origtree := IxLxI_kmn_n( TL(16, 2, 1, 1).withTags([ AVecReg(ppu_4x32f) ]),
          SIMD_ISA_Bases2( TL(8, 2, 1, 2).withTags([ AVecReg(ppu_4x32f) ]) ),
          SIMD_ISA_Bases1( TL(4, 2, 4, 1).withTags([ AVecReg(ppu_4x32f) ]) ) ),
      measured := 0,
      globalUnrolling := 8 ) ]);
SIMD_ISA_DB.hashAdd(TL(8, 2, 1, 1).withTags([ AVecReg(spu_4x32f) ]), [ rec(
      ruletree := SIMD_ISA_Bases1( TL(8, 2, 1, 1).withTags([ AVecReg(ppu_4x32f) ]) ),
      origtree := SIMD_ISA_Bases1( TL(8, 2, 1, 1).withTags([ AVecReg(ppu_4x32f) ]) ),
      measured := 0,
      globalUnrolling := 8 ) ]);
SIMD_ISA_DB.hashAdd(TL(4, 2, 1, 1).withTags([ AVecReg(spu_2x64f) ]), [ rec(
      ruletree := SIMD_ISA_Bases1( TL(4, 2, 1, 1).withTags([ AVecReg(spu_2x64f) ]) ),
      origtree := SIMD_ISA_Bases1( TL(4, 2, 1, 1).withTags([ AVecReg(spu_2x64f) ]) ),
      measured := 0,
      globalUnrolling := 8 ) ]);
SIMD_ISA_DB.hashAdd(TL(1, 1, 1, 2).withTags([ AVecReg(spu_2x64f) ]), [ rec(
      ruletree := IxLxI_vtensor( TL(1, 1, 1, 2).withTags([ AVecReg(spu_2x64f) ]) ),
      origtree := IxLxI_vtensor( TL(1, 1, 1, 2).withTags([ AVecReg(spu_2x64f) ]) ),
      measured := 0,
      globalUnrolling := 8 ) ]);
SIMD_ISA_DB.hashAdd(TL(4, 2, 1, 2).withTags([ AVecReg(spu_4x32f) ]), [ rec(
      ruletree := SIMD_ISA_Bases1( TL(4, 2, 1, 2).withTags([ AVecReg(ppu_4x32f) ]) ),
      origtree := SIMD_ISA_Bases1( TL(4, 2, 1, 2).withTags([ AVecReg(ppu_4x32f) ]) ),
      measured := 0,
      globalUnrolling := 8 ) ]);
SIMD_ISA_DB.hashAdd(TL(4, 2, 2, 1).withTags([ AVecReg(spu_4x32f) ]), [ rec(
      ruletree := SIMD_ISA_Bases1( TL(4, 2, 2, 1).withTags([ AVecReg(ppu_4x32f) ]) ),
      origtree := SIMD_ISA_Bases1( TL(4, 2, 2, 1).withTags([ AVecReg(ppu_4x32f) ]) ),
      measured := 0,
      globalUnrolling := 8 ) ]);
SIMD_ISA_DB.hashAdd(TL(4, 2, 1, 4).withTags([ AVecReg(spu_4x32f) ]), [ rec(
      ruletree := IxLxI_vtensor( TL(4, 2, 1, 4).withTags([ AVecReg(ppu_4x32f) ]) ),
      origtree := IxLxI_vtensor( TL(4, 2, 1, 4).withTags([ AVecReg(ppu_4x32f) ]) ),
      measured := 0,
      globalUnrolling := 8 ) ]);
SIMD_ISA_DB.hashAdd(TL(4, 2, 4, 1).withTags([ AVecReg(spu_4x32f) ]), [ rec(
      ruletree := SIMD_ISA_Bases1( TL(4, 2, 4, 1).withTags([ AVecReg(ppu_4x32f) ]) ),
      origtree := SIMD_ISA_Bases1( TL(4, 2, 4, 1).withTags([ AVecReg(ppu_4x32f) ]) ),
      measured := 0,
      globalUnrolling := 8 ) ]);
SIMD_ISA_DB.hashAdd(TL(16, 8, 1, 1).withTags([ AVecReg(spu_4x32f) ]), [ rec(
      ruletree := IxLxI_kmn_km( TL(16, 8, 1, 1).withTags([ AVecReg(ppu_4x32f) ]),
          SIMD_ISA_Bases1( TL(4, 2, 4, 1).withTags([ AVecReg(ppu_4x32f) ]) ),
          SIMD_ISA_Bases2( TL(8, 4, 1, 2).withTags([ AVecReg(ppu_4x32f) ]) ) ),
      origtree := IxLxI_kmn_km( TL(16, 8, 1, 1).withTags([ AVecReg(ppu_4x32f) ]),
          SIMD_ISA_Bases1( TL(4, 2, 4, 1).withTags([ AVecReg(ppu_4x32f) ]) ),
          SIMD_ISA_Bases2( TL(8, 4, 1, 2).withTags([ AVecReg(ppu_4x32f) ]) ) ),
      measured := 0,
      globalUnrolling := 10000 ) ]);
SIMD_ISA_DB.hashAdd(TL(8, 4, 1, 1).withTags([ AVecReg(spu_4x32f) ]), [ rec(
      ruletree := SIMD_ISA_Bases1( TL(8, 4, 1, 1).withTags([ AVecReg(ppu_4x32f) ]) ),
      origtree := SIMD_ISA_Bases1( TL(8, 4, 1, 1).withTags([ AVecReg(ppu_4x32f) ]) ),
      measured := 0,
      globalUnrolling := 8 ) ]);
SIMD_ISA_DB.hashAdd(TL(8, 4, 1, 2).withTags([ AVecReg(spu_4x32f) ]), [ rec(
      ruletree := SIMD_ISA_Bases2( TL(8, 4, 1, 2).withTags([ AVecReg(ppu_4x32f) ]) ),
      origtree := SIMD_ISA_Bases2( TL(8, 4, 1, 2).withTags([ AVecReg(ppu_4x32f) ]) ),
      measured := 0,
      globalUnrolling := 8 ) ]);
SIMD_ISA_DB.hashAdd(TL(8, 4, 2, 1).withTags([ AVecReg(spu_4x32f) ]), [ rec(
      ruletree := SIMD_ISA_Bases1( TL(8, 4, 2, 1).withTags([ AVecReg(ppu_4x32f) ]) ),
      origtree := SIMD_ISA_Bases1( TL(8, 4, 2, 1).withTags([ AVecReg(ppu_4x32f) ]) ),
      measured := 0,
      globalUnrolling := 8 ) ]);
SIMD_ISA_DB.hashAdd(TL(8, 2, 1, 2).withTags([ AVecReg(spu_4x32f) ]), [ rec(
      ruletree := SIMD_ISA_Bases2( TL(8, 2, 1, 2).withTags([ AVecReg(ppu_4x32f) ]) ),
      origtree := SIMD_ISA_Bases2( TL(8, 2, 1, 2).withTags([ AVecReg(ppu_4x32f) ]) ),
      measured := 0,
      globalUnrolling := 8 ) ]);
SIMD_ISA_DB.hashAdd(TL(8, 2, 2, 1).withTags([ AVecReg(spu_4x32f) ]), [ rec(
      ruletree := SIMD_ISA_Bases1( TL(8, 2, 2, 1).withTags([ AVecReg(ppu_4x32f) ]) ),
      origtree := SIMD_ISA_Bases1( TL(8, 2, 2, 1).withTags([ AVecReg(ppu_4x32f) ]) ),
      measured := 0,
      globalUnrolling := 8 ) ]);
SIMD_ISA_DB.hashAdd(TL(16, 4, 1, 1).withTags([ AVecReg(spu_4x32f) ]), [ rec(
      ruletree := IxLxI_kmn_n( TL(16, 4, 1, 1).withTags([ AVecReg(ppu_4x32f) ]),
          SIMD_ISA_Bases2( TL(8, 4, 1, 2).withTags([ AVecReg(ppu_4x32f) ]) ),
          SIMD_ISA_Bases1( TL(8, 4, 2, 1).withTags([ AVecReg(ppu_4x32f) ]) ) ),
      origtree := IxLxI_kmn_n( TL(16, 4, 1, 1).withTags([ AVecReg(ppu_4x32f) ]),
          SIMD_ISA_Bases2( TL(8, 4, 1, 2).withTags([ AVecReg(ppu_4x32f) ]) ),
          SIMD_ISA_Bases1( TL(8, 4, 2, 1).withTags([ AVecReg(ppu_4x32f) ]) ) ),
      measured := 0,
      globalUnrolling := 8 ) ]);
SIMD_ISA_DB.hashAdd(TL(4, 2, 2, 2).withTags([ AVecReg(spu_4x32f) ]), [ rec(
      ruletree := SIMD_ISA_Bases1( TL(4, 2, 2, 2).withTags([ AVecReg(ppu_4x32f) ]) ),
      origtree := SIMD_ISA_Bases1( TL(4, 2, 2, 2).withTags([ AVecReg(ppu_4x32f) ]) ),
      measured := 0,
      globalUnrolling := 8 ) ]);
SIMD_ISA_DB.hashAdd(TL(16, 2, 1, 1).withTags([ AVecReg(spu_4x32f) ]), [ rec(
      ruletree := IxLxI_kmn_n( TL(16, 2, 1, 1).withTags([ AVecReg(ppu_4x32f) ]),
          SIMD_ISA_Bases2( TL(8, 2, 1, 2).withTags([ AVecReg(ppu_4x32f) ]) ),
          SIMD_ISA_Bases1( TL(4, 2, 4, 1).withTags([ AVecReg(ppu_4x32f) ]) ) ),
      origtree := IxLxI_kmn_n( TL(16, 2, 1, 1).withTags([ AVecReg(ppu_4x32f) ]),
          SIMD_ISA_Bases2( TL(8, 2, 1, 2).withTags([ AVecReg(ppu_4x32f) ]) ),
          SIMD_ISA_Bases1( TL(4, 2, 4, 1).withTags([ AVecReg(ppu_4x32f) ]) ) ),
      measured := 0,
      globalUnrolling := 8 ) ]);
SIMD_ISA_DB.hashAdd(TL(8, 2, 1, 1).withTags([ AVecReg(spu_4x32f) ]), [ rec(
      ruletree := SIMD_ISA_Bases1( TL(8, 2, 1, 1).withTags([ AVecReg(ppu_4x32f) ]) ),
      origtree := SIMD_ISA_Bases1( TL(8, 2, 1, 1).withTags([ AVecReg(ppu_4x32f) ]) ),
      measured := 0,
      globalUnrolling := 8 ) ]);


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

CoSynthesize := function(t, opts)
    local tmp, f;
    tmp := Copy(t);
    for f in opts.csStrategy do
        tmp := f(tmp, opts);
    od;
    return tmp;
end;

ExportCode := function(c)
    local vars;
    vars := FoldL(Collect(c, var), (a,b)->When(not b.id in List(a, i->i.id), Concat([b],a), a), []);
    Print("let(", DoForAll(vars, i->Print(i.id, " := var(\"", i.id, "\", ", i.t, "),\n")), c, ")");
end;

ExportCProg := function(c, opts)
    Print(opts.doc::
        "\n");
    PrintCode(opts.subName, c, opts);
end;

ExportInclude := function(c, opts)
    local parameters, id, unparser, oopts, macro;

    unparser := Copy(opts.unparser);
    opts.unparser.opts:=opts;
    macro := "__"::StringToUpper(opts.filename)::"_H__";

    Print(opts.doc::
        "\n"::
        "#ifndef "::macro::"\n"::
        "#define "::macro::"\n\n"::
        "#ifdef __cplusplus\n"::
        "extern \"C\" {\n"::
        "#endif\n");

    parameters:=Flat(c.params);
    id := opts.subName;
    Print("\n",
        When(IsBound(c.inline) and c.inline, "inline ",""),
        opts.funcModifier, opts.unparser.declare(c.ret, var(id, c.ret), 0, 1), "(",
            DoForAllButLast(parameters, p->Print(unparser.declare(p.t, p,0,1), ", ")),
            When(Length(parameters)>0, unparser.declare(Last(parameters).t, Last(parameters),0,1), ""), ");\n");

    Print("\n#ifdef __cplusplus\n"::
        "}\n"::
        "#endif\n\n"::
        "#endif\n");
    opts.unparser := unparser;
end;

max_stars_repo_name: spiral-software/spiral-package-hcol

max_issues_repo_name: spiral-software/spiral-package-hcol

max_forks_repo_name: spiral-software/spiral-package-hcol


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

# Cross-compatibility
#

#######################################################################################
#   tSPL DFT rules
NewRulesFor(DFT, rec(
    #F DFT_CT: 1965
    #F   General Cooley-Tukey Rule
    #F   DFT_n = (DFT_n/d tensor I_d) * diag * (I_n/d tensor F_d) * perm
    #F
    #F Cooley/Tukey:
    #F   An Algorithm for the Machine Calculation of Complex Fourier Series.
    #F   Mathematics of Computation, Vol. 19, 1965, pp. 297--301.
    #F
    DFT_tSPL_CT := rec(
    info          := "tSPL DFT(mn,k) -> DFT(m, k%m), DFT(n, k%n)",

    maxSize       := false,
    filter := e->true,

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

    children      := (self, nt) >> Map2(Filtered(DivisorPairs(nt.params[1]), self.filter), (m,n) -> [
        TCompose([
            TGrp(TCompose([
                TTensorI(DFT(m, nt.params[2] mod m), n, AVec, AVec),
                TTwiddle(m*n, n, nt.params[2])
            ])),
            TGrp(TTensorI(DFT(n, nt.params[2] mod n), m, APar, AVec))
        ]).withTags(nt.getTags())
    ]),

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

#D    isApplicable := (self,P) >> P[1] > 2 and
#D        (self.maxSize=false or P[1] <= self.maxSize) and not IsPrime(P[1]) and PHasTags(self.nonTerminal, P),
#D    allChildren  := P -> Map2(DivisorPairs(P[1]),
#D        (m,n) -> [ TCompose([TGrp(TCompose([TTensorI(DFT(m, P[2] mod m), n, AVec, AVec), TDiag(fPrecompute(Tw1(m*n, n, P[2])))])),
#D                   TGrp(TTensorI(DFT(n, P[2] mod n), m, APar, AVec))], P[3]) ]),
#D    rule := (P,C,nt) -> C[1],
    switch := false
    )
));

NewRulesFor(DFT, rec(
    DFT_TTensorI_CT := rec(
        switch := false,
        maxSize := false,
        minSize := false,

        applicable := (self, t) >> let(
            n := Rows(t),
            n > 2
            #and nt.hasTags() --> fixed: dpickem 08/03/09
            and t.hasTags()
            and (self.maxSize=false or n <= self.maxSize)
            and (self.minSize=false or n >= self.minSize)
            and not IsPrime(n)
        ),

        children := t -> let(
            tags := t.params[3],
            Map2(DivisorPairs(Rows(t)), (m,n) -> [
                TTensorI_GT(GT(DFT(m, t.params[2] mod m), XChain([0, 1]), XChain([0, 1]), [n]).withTags(tags)),
                TTensorI_GT(GT(DFT(n, t.params[2] mod n), XChain([0, 1]), XChain([1, 0]), [m]).withTags(tags))
            ])
        ),

        apply := (t, C, Nonterms) ->
            Grp(C[1] * Diag(fPrecompute(Tw1(Rows(t), Rows(Nonterms[2].params[1]), t.params[2])))) * C[2]
    ),

    #F DFT_Rader: Rader's Algorithm for Prime size DFT
    #F
    #F   DFT(p) -> P(g)' * (1 dirsum DFT(p-1)) * Tp * (1 dirsum DFT(p-1)) * P(g)",
    #F   P(g) = perm
    #F   Tp   = [[1,1],[1,-1/(p-1)]] dirsum diag
    #F
    DFT_tSPL_Rader := rec(
        minSize := 3,
        useSymmetricAlgorithm := true,
        avoidSizes := [],
        switch := false,

        applicable := (self, t) >> let(
            n := Rows(t),
            n >= self.minSize
            and IsPrime(n)
            and (not n in self.avoidSizes)
            and t.hasTags()
        ),

        children := (self, t) >> let(
            N := Rows(t),
            tags := t.getTags(),
            When(self.useSymmetricAlgorithm,
                [[
                    TDirectSum(I(1), DFT(N-1, -1)).withTags(tags),
                    TRaderMid(N, t.params[2], PrimitiveRootMod(N)).withTags(tags)
                ]],
                [[
                    TCompose([
                        TDirectSum(I(1), DFT(N-1, -1)),
                        TRaderMid(N, t.params[2], PrimitiveRootMod(N)),
                        TDirectSum(I(1), DFT(N-1, -1))
                    ]).withTags(tags)
                ]]
            )
        ),

    apply := (self, t, C, Nonterms) >> let(
        N := Rows(t), k := t.params[2], root := PrimitiveRootMod(N),
        When(self.useSymmetricAlgorithm,
            RR(N, 1, root).transpose()
            * C[1].transpose() * C[2] * C[1]
            * RR(N, 1, root),
            RR(N, 1, root).transpose()
            * C[1]
            * RR(N, 1, root)
        ))
    ),

    #F DFT_Bluestein : Convert FFT to Toeplitz matrix, then embed in a larger Circulant
    #F
    #F    DFT_n -> diag * Toeplitz * diag
    #F
    # DFT on inputs, some of which are 0:
    #
    DFT_tSPL_Bluestein := rec(
        minRoundup := 8,
        customFilter := True,
        goodFactors := [2,3,5],
        applicableSizes := IsPrime,
        circSizes := meth(self,N)
            local low, high, cands;
            low := 2 * N + 1;
            high := 4 * 2^Log2Int(N);
            cands := When(high < self.minRoundup, [self.minRoundup], [low..high]);
            cands := Filtered(cands, self.customFilter);
            cands := Filtered(cands, i->IsEvenInt(i) and IsSubset(Set(self.goodFactors), Set(Factors(i))));
            return cands;
        end,

        toeplitz := (N, k) ->
            List([1..2*N-1], x -> ComplexW(2*N, -(N-x)^2 * k)),

        diag := (N, k) ->
            List([0..N-1], x -> ComplexW(2*N, x^2 * k)),

        circulantToep := (toep,size) -> let(l:=Length(toep),
            Concatenation(
                List([1..(l+1)/2],   x->toep[(l+1)/2+x-1]),
                List([1..size-l],    x->0),
                List([1..(l-1)/2],   x->toep[x]))),

        switch := false,
        applicable := (self, t) >> self.applicableSizes(Rows(t)) and t.hasTags(),

        children := meth(self, t)
            local P, N, k, circsize, dvar, diag, toep, circ, diagonalized_circ, tags, kids;

            tags := t.getTags();
            P := t.params;
            N := P[1];
            k := P[2];
            diag := FData(self.diag(N,k));
            toep := self.toeplitz(N,k);
            kids := [];
            for circsize in self.circSizes(N) do
                dvar := Dat(TArray(TComplex, circsize));
                circ := self.circulantToep(toep, circsize);
                diagonalized_circ := FData(1/circsize * ComplexFFT(circ));
                Add(kids, [
                    TCompose([
                        TRDiag(circsize/2, N, diag).transpose(),
                        PrunedDFT(circsize, -1, circsize/2, [0]).transpose(),
                        TDiag(diagonalized_circ),
                        PrunedDFT(circsize, 1, circsize/2, [0]),
                        TRDiag(circsize/2, N, diag)
                    ]).withTags(tags)
                ]);
            od;

            return kids;
        end,

        apply := (t, C, Nonterms) -> C[1]
    ),

    #F DFT_GoodThomas : Prime Factor FFT
    #F
    #F     DFT_n*k -> perm * (DFT_n_a tensor DFT_k_b) * perm
    #F     when gcd(n,k) = 1
    #F
    DFT_tSPL_GoodThomas := rec(
        applicable     := (t) ->
            t.params[1] > 2
            and DivisorPairsRP(t.params[1]) <> []
            and t.hasTags(),

        children := (self, t) >> let(
            tags := t.getTags(),
            mn := t.params[1],
            k := t.params[2],
            Map2(DivisorPairsRP(mn), (m,n) -> [
                TTensor(DFT(m, k*n mod m), DFT(n, k*m mod n)).withTags(tags)
            ])
        ),

        apply := (t, C, Nonterms) -> let(
            P := t.params,
            r := Rows(Nonterms[1].params[1]),
            s := Rows(Nonterms[1].params[2]),
            alpha := 1 / s mod r,
            beta  := 1 / r mod s,
            # CRT should not be pulled into sums before DelayedDirectSums are terminated,
            # otherwise PFA inside Rader for n =vq breaks as VectRaderDiag rules can't match
            DelayedPrm(CRT(r,s,1,1)).transpose() * C[1] * DelayedPrm(CRT(r,s,1,1))
        )
    )
));

#--------------------------------------------------------------------------------------

NewRulesFor(TTensorI, rec(

    # Vector recursion for the DFT
    DFT_vecrec := rec(
        forTransposition := false,
        applicable := nt ->
            ObjId(nt.params[1]) = DFT
            and IsParVec(nt.params)
            and let(p1 := nt.params[1].params[1],
                p1 > 2 and not IsPrime(p1)
            )
            and nt.hasTags(),

        children := nt -> let(
            f := nt.params[1],
            k := nt.params[2],
            Map2(DivisorPairs(f.params[1]), (m,n) -> let(
                r := When(
                    nt.hasTags() and nt.firstTagIs(spiral.paradigms.vector.AVecReg),
                    VWrapTRC(nt.firstTag().v),
                    VWrapId
                ),
                [
                    TGrp(TCompose([
                        TTensorI(DFT(m, f.params[2] mod m), n, AVec, AVec),
                        TDiag(fPrecompute(Tw1(m*n, n, f.params[2])))
                    ])).withTags(nt.getTags()), <# strange #>
                    TTensorI(
                        TTensorI(DFT(n, f.params[2] mod n), m, APar, AVec),
                        k, APar, AVec
                    ).withTags(nt.getTags()) <# strange #>
                ]
            ))
        ),

        apply := (nt, c, cnt) -> Tensor(I(nt.params[2]), c[1]) * c[2],

        switch := false

#D    isApplicable := (self,P) >> ObjId(P[1])=DFT and P[3].isPar and P[4].isVec and
#D                                    let(P1:=P[1].params[1], P1 > 2 and not IsPrime(P1)) and PHasTags(self.nonTerminal, P),

#D        allChildren  := P -> let(k:=P[2],
#D        Map2(DivisorPairs(P[1].params[1]),
#D        (m,n) -> let( r:=When(Length(P[5])>0 and P[5][1].isReg, VWrapTRC(P[5][1].v), VWrapId),
#D            [ TGrp(TCompose([TTensorI(DFT(m, P[1].params[2] mod m), n, AVec, AVec),
#D                     TDiag(fPrecompute(Tw1(m*n, n, P[1].params[2])))]), P[5]).setWrap(r),

#D              TTensorI(TTensorI(DFT(n, P[1].params[2] mod n), m, APar, AVec),
#D                  k, APar, AVec, P[5]).setWrap(r) ]))),

#D    rule := (P,C,nt) -> Tensor(I(P[2]), C[1]) * C[2],
    ),

    # Transposed vector recursion for the DFT
    DFT_vecrec_T := rec(
        forTransposition := false,

        applicable := nt ->
            ObjId(nt.params[1]) = DFT
            and IsVecPar(nt.params)
            and let(
                p1 := nt.params[1].params[1],
                p1 > 2
                and IsPrime(p1)
            )
            and nt.hasTags(),

        children := nt -> let(
            k := nt.params[2],
            Map2(DivisorPairs(nt.params[1].params[1], (m,n) -> let(
                r := When(nt.hasTags() and nt.isTag(1, spiral.paradigms.vector.AVecReg),
                    VWrapTRC(nt.firstTag().v),
                    VWrapId
                ),
                [
                    TTensorI(
                        TTensorI(DFT(n, P[1].params[2] mod n), m, AVec, APar),
                        k, AVec, APar
                    ).withTags(nt.getTags()).setWrap(r), <# strange #>
                    TGrp(TCompose([
                        TDiag(fPrecompute(Tw1(m*n, n, P[1].params[2]))),
                        TTensorI(DFT(m, P[1].params[2] mod m), n, AVec, AVec)
                    ])).withTags(nt.getTags()).setWrap(r) <# strange #>
                ]
            )))
        ),

        apply := (nt, c, cnt) -> c[1] * Tensor(I(nt.params[2]), c[2]),

        switch := false

#D        isApplicable := (self,P) >> ObjId(P[1])=DFT and P[3].isVec and P[4].isPar and
#D                                    let(P1:=P[1].params[1], P1 > 2 and not IsPrime(P1)) and PHasTags(self.nonTerminal, P),
#D        allChildren  := P -> let(k:=P[2],
#D            Map2(DivisorPairs(P[1].params[1]),
#D        (m,n) -> let(r:=When(Length(P[5])>0 and P[5][1].isReg, VWrapTRC(P[5][1].v), VWrapId), [
#D            TTensorI(TTensorI(DFT(n, P[1].params[2] mod n), m, AVec, APar),
#D                                        k, AVec, APar, P[5]).setWrap(r),
#D
#D                    TGrp(TCompose([TDiag(fPrecompute(Tw1(m*n, n, P[1].params[2]))),
#D                                   TTensorI(DFT(m, P[1].params[2] mod m), n, AVec, AVec)
#D                  ]), P[5]).setWrap(r)
#D                            ]))),
#D
#D    rule := (P,C) -> C[1] * Tensor(I(P[2]), C[2]),

    )
));

#######################################################################################################
#   tSPL rule
NewRulesFor(MDDFT, rec(
    MDDFT_tSPL_RowCol := rec(
        info := "tSPL MDDFT_n -> MDDFT_n/d, MDDFT_d",

        applicable := (self, t) >> Length(t.params[1]) > 1,
        freedoms := t -> [ [1..Length(t.params[1])-1] ],

        child := (t, fr) -> let(
            newdims := SplitAt(t.params[1], fr[1]),
            rot := t.params[2],
            [ TTensor(
                MDDFT(newdims[1], rot),
                MDDFT(newdims[2], rot)
            ).withTags(t.getTags())]
        ),

        apply := (t, C, Nonterms) -> C[1],
        switch := false
    ),
));

AllRadices := function(n)
    local divisors, radices;
    divisors := DropLast(Drop(DivisorsInt(n), 1), 1);
    radices := Filtered(divisors, i -> n=(i^LogInt(n, i)));
    return radices;
end;

##########################################################################################
#   tSPL Pease DFT rule
NewRulesFor(DFT, rec(
    DFT_tSPL_Pease   := rec (
        info             := "Pease tSPL DFT_(k) -> (\Prod(L(I tensor F)Tc))*DR(k,r)",
#        forTransposition := false,
        forTransposition := true,
        maxRadix         := 32,
        minRadix         := 4,

        applicable := (self, nt) >>
            nt.hasTags()
            and nt.params[2] = 1
            and let(
                r := AllRadices(nt.params[1]),
                Length(r) > 0
                and ForAny(r, i ->
                    i >= self.minRadix
                    and i <= self.maxRadix
                    #and nt.firstTag().legal_kernel(i)
                )
            ),

        children := (self, nt) >> let(
            N := nt.params[1],
            m := nt.params[2],
            radices := Filtered(AllRadices(N), i ->
                i >= self.minRadix and i <= self.maxRadix
                #and nt.firstTag().legal_kernel(i)
            ),
            j := var("j"),

            List(radices, rdx -> [
                TCompose([
                    TICompose(j, LogInt(N, rdx),
                        TCompose([
                            TTensorI(DFT(rdx, m), N/rdx, AVec, APar),
                            TDiag(fPrecompute(TC(N, rdx, j, m)))
                        ])
                    ),
                    TDR(N, rdx)
                ]).withTags(nt.getTags()),
            ])
        ),

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

        switch := false

#D    isApplicable     := (self, P) >> PHasTags(self.nonTerminal, P) and P[2] = 1 and
#D                                     let(r := AllRadices(P[1]),
#D                                        Length(r) > 0 and
#D                                        ForAny(r, i-> i >= self.minRadix and i <= self.maxRadix and PGetFirstTag(self.nonTerminal, P).legal_kernel(i))
#D                                     ),
#D
#D    allChildren      := (self, P) >> let(
#D                            N := P[1],
#D                            m := P[2],
#D                            radices := Filtered(AllRadices(N), i-> i >= self.minRadix and i <= self.maxRadix and PGetFirstTag(self.nonTerminal, P).legal_kernel(i)),
#D                            j := var("j"),
#D
#D                            List(radices, rdx ->
#D                                [
#D                                    AddTag(
#D                                        TCompose([
#D                                            TICompose(j, LogInt(N, rdx),
#D                                                TCompose([
#D                                                    TTensorI(DFT(rdx, m), N/rdx, AVec, APar),
#D                                                    TDiag(fPrecompute(TC(N, rdx, j, m)))
#D                                                ])),
#D                                            TDR(N, rdx)
#D                                        ]),
#D                                    PGetTags(self.nonTerminal, P))
#D                                ]
#D                            )),
#D
#D    rule := (P, C) -> C[1],
    )
));

##########################################################################################
#   tSPL Stockham DFT rule
NewRulesFor(DFT, rec(
    DFT_tSPL_Stockham   := rec (
        info             := "Stockham tSPL DFT_(k) -> (\Prod(DFT tensor I)*Diag*(L tensor I))",
        forTransposition := false,
        maxRadix         := 2,
        minRadix         := 2,

        #one dimensional transform --> nt.params[2] = 1
        applicable := (self, nt) >>
            nt.hasTags()
            and nt.params[2] = 1
            and let(
                r := AllRadices(nt.params[1]),
                Length(r) > 0
                and ForAny(r, i ->
                    i >= self.minRadix
                    and i <= self.maxRadix
                    #and nt.firstTag().legal_kernel(i)
                )
            ),

        children := (self, nt) >> let(
            N := nt.params[1],
            m := nt.params[2],
            radices := Filtered(AllRadices(N), i ->
                i >= self.minRadix and i <= self.maxRadix and (When(IsBound(nt.firstTag().legal_kernel), nt.firstTag().legal_kernel(i), true))
            ),
#need the tag for TICompose to be parallelized, but then the DFT(r) in TTensorI will not be broken down
#if i use no tag for TICompose, no parallelization will take place ...
#solution: drop tags in GT_Par's children function --> rule GT_Par_drop implemented in paradigms/gpu/breakdown.gi
            List(radices, rdx -> let(j := var.fresh("j", TInt, LogInt(N, rdx)), tags := Drop(t.getTags(), 1),
                [TICompose(j, LogInt(N, rdx),
                        TCompose([
                            TTensorI(DFT(rdx, m), N/rdx, AVec, AVec),
                            #Diag(fPrecompute(diagTensor(Stockham_radix.gen(rdx, LogInt(N, rdx)-j-1), fConst(TComplex, rdx^j, 1)))),
                            #Tensor(Diag(fPrecompute(Stockham_radix.gen(rdx, LogInt(N, rdx)-j-1))), I(rdx^j)),
                            #TDiag(fPrecompute(diagTensor(TC5(rdx, LogInt(N, rdx)-j-1), fConst(TComplex, rdx^j, 1)))),
                            TTwiddle_Stockham(N, N/rdx, rdx, j),
                            Tensor(L(N/rdx^j, rdx), I( rdx^j))
                        ])
                    ).withTags(nt.getTags()),
                ])
            )
        ),

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

        switch := false
    )
));

NewRulesFor(DFT, rec(
    DFT_tSPL_Stockham_split   := rec (
        info             	:= "Stockham tSPL DFT_(k) -> (\Prod(DFT tensor I)*Diag*(L tensor I)), first n iterations split off",
        forTransposition 	:= false,
        minRadix         	:= 4,
        maxRadix         	:= 4,

        # one dimensional transform --> nt.params[2] = 1
        # first parameter of AGenericTag indicates how many iterations to split off
        applicable := (self, nt) >>
            nt.hasTags()
#            and nt.firstTag().kind() = AGenericTag
#            and IsPosInt(nt.firstTag().params[1])
            and nt.params[2] = 1
            and let(
                r := AllRadices(nt.params[1]),
                Length(r) > 0
                and ForAny(r, i ->
                    i >= self.minRadix
                    and i <= self.maxRadix
                    #and nt.firstTag() is GPU tag
                    #and nt.firstTag().legal_kernel(i)
                )
            ),

        children := (self, nt) >> let(
            N := nt.params[1],
            m := nt.params[2],
# number of iterations to peel off is n
#            n := nt.firstTag().params[1],
            n := 1,
            radices := Filtered(AllRadices(N), i ->
                i >= self.minRadix and i <= self.maxRadix and (When(IsBound(nt.firstTag().legal_kernel), nt.firstTag().legal_kernel(i), true))
            ),
            List(radices, rdx -> let(j := var.fresh("j", TInt, LogInt(N, rdx)-n), i := var.fresh("i", TInt, n),
            	# first TICompose deals with the iterations [n..LogInt(N,rdx)-1], these are the critical write-out stages

            	[TCompose([
                	TICompose(j, LogInt(N, rdx)-n,
                    	    TCompose([
                        	    TTensorI(DFT(rdx, m), N/rdx, AVec, AVec),
                            	#TDiag(fPrecompute(diagTensor(Stockham_radix.gen(rdx, LogInt(N, rdx)-j-1), fConst(TComplex, rdx^j, 1)))),
                            	#Diag(fPrecompute(diagTensor(TTwiddle_Stockham(rdx^(LogInt(N, rdx)-j), rdx^(LogInt(N, rdx)-j-1), 1).terminate(), fConst(TComplex, rdx^j, 1)))),
	                            TTwiddle_Stockham(N, N/rdx, rdx, j),
	                            TL(N/rdx^j, rdx, 1, rdx^j)
    	                    ])
        	            ).withTags(nt.getTags()),
            	     # second TICompose has to range from [0..n-1]
                	 TICompose(i, n,
	                        TCompose([
    	                        TTensorI(DFT(rdx, m), N/rdx, AVec, AVec),
        	                    #TDiag(fPrecompute(diagTensor(Stockham_radix.gen(rdx, n-i-1), fConst(TComplex, N/rdx^(n-i), 1)))),
                            	#Tensor(TTwiddle(rdx^(n-j), rdx^(n-j-1)), I(rdx^j)),
	                            #Diag(fPrecompute(diagTensor(TTwiddle_Stockham(rdx^(LogInt(N, rdx)-j), rdx^(LogInt(N, rdx)-j-1), 1).terminate(), fConst(TComplex, rdx^j, 1)))),
								TTwiddle_Stockham(N, N/rdx, rdx, j),
            	                TL(rdx^(n-i), rdx, 1, N/rdx^(n-i))
                	        ])
	                    ).withTags(nt.getTags()),
					])
                ])
            )
        ),

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

        switch := false
    )
));


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

BindGlobal("AG_SemigroupIterator",
function(arg)
  local next_iterator, is_done_iterator, shallow_copy,
        create_iter_rec, G, max_len;

  max_len := infinity;
  G := arg[1];
  if Length(arg) > 1 then
    max_len := arg[2];
    if not IsMonoid(G) and IsInt(max_len) then
      max_len := max_len - 1;
    fi;
  fi;

  create_iter_rec := function(G)
    local record, gens;

    record := rec(
      NextIterator := next_iterator,
      IsDoneIterator := is_done_iterator,
      ShallowCopy := shallow_copy,

      group := G,
      done := false,
      max_len := max_len,

      gens_ptr := [1, 1],

      ptr := 1,
    );

    if IsGroup(G) then
      # GAP is smart, it knows that a semigroup generated by group
      # generators of a finite group is in fact the whole group
      gens := GeneratorsOfGroup(G);
      record.gens := Difference(Concatenation(gens, List(gens, g->g^-1)), [One(G)]);
      record.elms := Concatenation([One(G)], record.gens);
      record.levels := [[1], [2..1+Length(record.gens)], []];
    elif IsMonoid(G) then
      record.gens := Difference(GeneratorsOfSemigroup(G), [One(G)]);
      record.elms := Concatenation([One(G)], record.gens);
      record.levels := [[1], [2..1+Length(record.gens)], []];
    else
      record.gens := ShallowCopy(GeneratorsOfSemigroup(G));
      record.elms := ShallowCopy(record.gens);
      record.levels := [[1..Length(record.gens)], []];
    fi;
    record.n_gens := Length(record.gens);

    if record.n_gens = 0 then
      record.is_list := true;
    elif IsMonoid(G) then
      if max_len = 0 then
        record.is_list := true;
        record.elms := [One(G)];
      else
        record.is_list := max_len = 1;
      fi;
    else
      if max_len = -1 then
        record.is_list := true;
        record.elms := [];
      else
        record.is_list := max_len = 0;
      fi;
    fi;

    return record;
  end;

  next_iterator := function(iter)
    local elm;

    if is_done_iterator(iter) then
      Error("<iter> is exhausted");
    fi;

    elm := iter!.elms[iter!.ptr];
    iter!.ptr := iter!.ptr + 1;

    return elm;
  end;

  is_done_iterator := function(iter)
    local gens_ptr, gens, elm, elm_list, levels, n_levels, added;

    if iter!.done then
      return true;
    fi;

    if iter!.ptr <= Length(iter!.elms) then
      return false;
    fi;

    if iter!.is_list then
      iter!.done := true;
      return true;
    fi;

    gens_ptr := iter!.gens_ptr;
    gens := iter!.gens;
    elm_list := iter!.elms;
    levels := iter!.levels;
    n_levels := Length(levels);
    added := false;

    while not added do
      if gens_ptr[1] > Length(levels[n_levels-1]) then
        if IsEmpty(levels[n_levels]) then
          iter!.done := true;
          SetIsFinite(iter!.group, true);
          SetSize(iter!.group, Length(elm_list));
          break;
        elif n_levels > iter!.max_len then
          iter!.done := true;
          break;
        fi;

        Add(levels, []);
        n_levels := n_levels + 1;
        gens_ptr[1] := 1;
      fi;

      elm := elm_list[levels[n_levels-1][gens_ptr[1]]] * gens[gens_ptr[2]];

      if IsGroup(iter!.group) then
        if not elm in elm_list{levels[n_levels-2]} and
           not elm in elm_list{levels[n_levels-1]} and
           not elm in elm_list{levels[n_levels]}
        then
          if elm in elm_list then
            Print("elm: ", elm, "\n");
            Print("elm_list: ", elm_list, "\n");
            Print("levels: ", levels, "\n");
            Print("n_levels: ", n_levels, "\n");
            Error("WTF");
          fi;
          Add(elm_list, elm);
          Add(levels[n_levels], Length(elm_list));
          added := true;
        fi;
      elif not elm in elm_list then
        Add(elm_list, elm);
        Add(levels[n_levels], Length(elm_list));
        added := true;
      fi;

      gens_ptr[2] := gens_ptr[2] + 1;
      if gens_ptr[2] > Length(gens) then
        gens_ptr[2] := 1;
        gens_ptr[1] := gens_ptr[1] + 1;
      fi;
    od;

    return iter!.done;
  end;

  shallow_copy := function(iter)
    return create_iter_rec(iter!.group);
  end;

  return IteratorByFunctions(create_iter_rec(G));
end);

################################################################################
##
#M  Iterator( <G>[, <max_len>] )
##
##  Provides a possibility to loop over elements of a group (semigroup, monoid)
##  generated by automata. If <max_len>  is given stops after enumerating all elements
##  of length up to <max_len>.
##  \beginexample
##  gap> Grigorchuk_Group := AutomatonGroup("a=(1,1)(1,2),b=(a,c),c=(a,d),d=(1,b)");
##  < a, b, c, d >
##  gap> iter := Iterator(Grigorchuk_Group,5);
##  <iterator>
##  gap> l := [];;
##  gap> for g in iter do
##  >      if Order(g) = 16 then Add(l,g); fi;
##  >    od;
##  gap> l;
##  [ b*a, a*b, d*a*c, c*a*d, d*a*c*a, d*a*b*a, c*a*d*a, b*a*d*a, a*d*a*c,
##    a*d*a*b, a*c*a*d, a*b*a*d, c*a*c*a*b, c*a*b*a*b, b*a*c*a*c, b*a*b*a*c,
##    a*d*a*c*a, a*c*a*d*a ]
##  \endexample
##

InstallMethod(Iterator, "for [IsTreeHomomorphismSemigroup]", [IsTreeHomomorphismSemigroup],
function(G)
  if IsGroup(G) and HasIsFinite(G) and IsFinite(G) then
    TryNextMethod();
  else
    return AG_SemigroupIterator(G, infinity);
  fi;
end);

InstallOtherMethod(Iterator, [IsTreeHomomorphismSemigroup, IsCyclotomic],
function(G, max_len)
  return AG_SemigroupIterator(G, max_len);
end);

###############################################################################
##
#M  TransformationSemigroupOnLevelOp (<G>, <k>)
##
InstallMethod(TransformationSemigroupOnLevelOp, "for [IsTreeHomomorphismSemigroup, IsPosInt]",
              [IsTreeHomomorphismSemigroup, IsPosInt],
function(G, k)
  local pgens, psemigroup, i;
  pgens := List(GeneratorsOfSemigroup(G), g -> TransformationOnLevel(g, k));

  # if there are Permutations and transformations convert all to transformations
  if Position(List(pgens,IsPerm),false) <> fail and Position(List(pgens,IsPerm),true) <> fail then
    for i in [1..Length(pgens)] do
      if IsPerm(pgens[i]) then pgens[i]:=AsTransformation(pgens[i]); fi;
    od;
  fi;

  psemigroup := SemigroupByGenerators(pgens);
  return psemigroup;
end);

###############################################################################
##
#M  MarkovOperator(<G>, <lev>, <weights>)
##
InstallMethod(MarkovOperator, "for [IsTreeHomomorphismSemigroup, IsCyclotomic, IsList]",
              [IsTreeHomomorphismSemigroup, IsCyclotomic, IsList],
function(G, n, weights)
  local gens,R,indet;
  gens := ShallowCopy(GeneratorsOfSemigroup(G));
  if Length(weights)<>Length(gens) then Error("The number of weights must match the number of generators"); fi;
  if IsString(weights[1]) then
    R:=PolynomialRing(Integers,weights);
    indet:=IndeterminatesOfPolynomialRing(R);
    return Sum(List([1..Length(gens)], x -> indet[x]*TransformationOnLevelAsMatrix(gens[x], n)));
  else
    return Sum(List([1..Length(gens)], x -> weights[x]*TransformationOnLevelAsMatrix(gens[x], n)));
  fi;
end);

###############################################################################
##
#M  MarkovOperator(<G>, <lev>)
##
InstallMethod(MarkovOperator, "for [IsTreeHomomorphismSemigroup, IsCyclotomic]",
              [IsTreeHomomorphismSemigroup, IsCyclotomic],
function(G, n)
  return MarkovOperator(G,n,List([1..Length(GeneratorsOfSemigroup(G))],x->1/(Length(GeneratorsOfSemigroup(G)))));
end);

###############################################################################
##
#M  Section(<G>, <v>)
##
InstallMethod(Section, "for [IsTreeHomomorphismSemigroup, IsList]",
              [IsTreeHomomorphismSemigroup, IsList],
function(G, v)
  local gens, g, sec_gens;
  gens:=GeneratorsOfSemigroup(G);
  for g in gens do
    if v^g<>v then
      Error("Section([IsTreeHomomorphismSemigroup, IsList]): <G> does not fix <v>");
      return fail;
    fi;
  od;

  if not IsAutomSemigroup(G) then
    return Semigroup(List(gens,x->Section(x,v)));
  else
    sec_gens:=[];
    for g in List(gens,x->Section(x,v)) do
      if not IsOne(g) and not g in sec_gens then
        Add(sec_gens, g);
      fi;
    od;
    if Length(sec_gens)>0 then
      return Semigroup(sec_gens);
    else
      return Semigroup([One(FamilyObj(Section(gens[1],v)))]);
    fi;
  fi;
end);

###############################################################################
##
#M  Section(<G>, <n>)
##
InstallMethod(Section, "for [IsTreeHomomorphismSemigroup, IsPosInt]",
              [IsTreeHomomorphismSemigroup, IsPosInt],
function(G, n)
  return Section(G,[n]);
end);

#E

max_stars_repo_name: gap-packages/automgrp

max_issues_repo_name: gap-packages/automgrp

max_forks_repo_name: gap-packages/automgrp


text:
QuadraticRoots := function(a, b, c)
  local d;
  d := Sqrt(b*b - 4*a*c);
  return [ (-b+d)/(2*a), (-b-d)/(2*a) ];
end;

# Hint : E(12) is a 12th primitive root of 1
QuadraticRoots(2, 2, -1);
# [ 1/2*E(12)^4-1/2*E(12)^7+1/2*E(12)^8+1/2*E(12)^11,
#   1/2*E(12)^4+1/2*E(12)^7+1/2*E(12)^8-1/2*E(12)^11 ]

# This works also with floating-point numbers
QuadraticRoots(2.0, 2.0, -1.0);
# [ 0.366025, -1.36603 ]


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

# Cell Unparser
# Author: schellap

@TInt := @.cond(x->x.t=TInt);
@TRealOrInt := @.cond(x->x.t=TInt or x.t=TReal);
@TReal := @.cond(x->x.t=TReal);
@TVect := @.cond(x->ObjId(x.t)=TVect);

Class(altivecUnparser, CMacroUnparserProg, rec(
    # str = format string with $n for arguments ($1 = first argument)
    # list of arguments , arguments are printed using this unparser
    printf := (self, str, args) >> ApplyFunc(PrintEvalF, Concatenation([str],
        List(args, a->(When(IsFunc(a), a, ()->self(a,0,0)))))),

    # -----------------------------
    # ISA independent constructs
    # -----------------------------
    nth :=  (self, o, i, is) >> self.printf("$1[$2]", [o.loc, o.idx]),
    fdiv := (self, o, i, is) >> self.printf("((($3)$1) / $2)", Concatenation(o.args, [self.opts.vector.isa.ctype])),
    fdiv := (self, o, i, is) >> self.printf("($1 / $2)", o.args),
    idiv := (self, o, i, is) >> self.printf("($1 / $2)", o.args),
    imod := (self, o, i, is) >> self.printf("($1 % $2)", o.args),

    # This is the type used for declarations of vector variables
    ctype := (t, isa) -> Cond(
        #t in [TDouble, TVect(TDouble, 1)], "float",
        t = TReal,              isa.ctype, #evals to "float" or "double"
        t = TVect(TReal, 2),    isa.vtype, #evals to "vector double"
        t = TVect(TReal, 4),    isa.vtype,
        t = TVect(TReal, 8),    isa.vtype,
        Error("Unparser doesn't know how to declare type: ", t)
    ),

    Value := (self, o, i, is) >> Cond(
        o.t = TReal,
           let(v := When(IsCyc(o.v), ReComplex(Complex(o.v)), Double(o.v)), When(v<0, Print("(", v, ")"), Print(v))),
        o.t = TInt,
           When(o.v < 0, Print("(", o.v, ")"), Print(o.v)),
        ObjId(o.t)=TVect,
           Print("((", self.ctype(o.t, self.opts.vector.isa), "){", self.infix(o.v, ", "), "})"),
        #NOTE: Adding this for parallel cell, but this shouldn't be needed, right?
       IsArray(o.t),
           Print("{", self.infix(o.v, ", "), "}")
    ),

    vdup         := (self, o, i, is) >> self.printf("vec_splat( ((vector float){$1, $1, $1, $1}), 0)", [o.args[1]]),

    # Declarations
    TVect := (self, t, vars, i, is) >> let(ctype := self.ctype(t, self.opts.vector.isa),
              Print(ctype, " ", self.infix(vars, ", "))),

    TReal := ~.TVect,

    TInt := (self, t, vars, i, is) >> Print("int ", self.infix(vars, ", ")),

    # Arithmetic
    # --------------------------------
    # -- mul --
    mul := (self, o, i, is) >> CondPat(o,
    [mul, @TReal, @TVect],
        Cond(self.opts.vector.isa.v = 2,
            self.printf("vec_madd(((vector $1){$2,$2}), $3, (vector $1)(0))",       [self.ctype(o.args[1].t, self.opts.vector.isa), o.args[1], o.args[2]]),
        self.opts.vector.isa.v = 4,
            self.printf("vec_madd(((vector $1){$2,$2,$2,$2}), $3, (vector $1)(0))", [self.ctype(o.args[1].t, self.opts.vector.isa), o.args[1], o.args[2]]),
        Error("Don't know how to unparse vector arch of length: ", self.opts.vector.isa.v)
        ),
    [mul, @TVect,   @TVect],
        self.printf("vec_madd($1, $2, (vector float)(0) )",              [o.args[1], o.args[2]]),
    [mul, @TRealOrInt, @TRealOrInt],
        self.printf("($1 * $2)", o.args),
    Error("Don't know how to unparse <o>. Unrecognized type combination")),

    # -- add --
    add := (self, o, i, is) >> When(Length(o.args) > 2,
       self(_computeExpType(add(o.args[1], _computeExpType(ApplyFunc(add, Drop(o.args, 1))))), i, is),
       CondPat(o,
           [add, @TVect,   @TVect],
              self.printf("vec_add($1, $2)", [o.args[1], o.args[2]]),
           #[add, @TRealOrInt, @TRealOrInt],
           [add, @, @],
              self.printf("($1 + $2)", o.args),
           Error("Don't know how to unparse <o>. Unrecognized type combination")
       )
    ),

    # -- sub --
    sub := (self, o, i, is) >> When(Length(o.args) > 2,
       self(_computeExpType(sub(o.args[1], _computeExpType(ApplyFunc(sub, Drop(o.args, 1))))), i, is),
       CondPat(o,
           [sub, @TVect,   @TVect],
              self.printf("vec_sub($1, $2)", [o.args[1], o.args[2]]),
           [sub, @TRealOrInt, @TRealOrInt],
              self.printf("($1 + $2)", o.args),
           Error("Don't know how to unparse <o>. Unrecognized type combination")
       )
    ),

    # -- neg --
    neg := (self, o, i, is) >> CondPat(o,
       [neg, @TVect],
           self(o.t.value(Replicate(o.t.size,0)) - o.args[1], i, is),
       self.printf("(-$1)", o.args)
    ),

    vpack := (self, o, i, is) >> Print(self.opts.vector.isa.vconstv, "{", self.infix(o.args, ", "), "}"),

    vunpackhi_4x32f_av := (self, o, i, is) >> self.printf("vec_mergeh($1, $2)", [o.args[1], o.args[2]]),
    vunpacklo_4x32f_av := (self, o, i, is) >> self.printf("vec_mergel($1, $2)", [o.args[1], o.args[2]]),

   vparam_av := (self, o, i, is) >>
      Print("((vector unsigned char){", PrintCS(prep_perm_string(o.p)), "})"),

   vzero_4x32f := (self, o, i, is) >> Print("((", self.opts.vector.isa.vtype ,"){", PrintCS(List([1..self.opts.vector.isa.v], i->0)), "})"),

   # ----------------------------------------------------------------------------------
   # ISA specific : altivec_4x32f
   # ----------------------------------------------------------------------------------
   vperm_4x32f := (self, o, i, is) >>
     self.printf("vec_perm($1, $2, $3)", o.args),

   vuperm_4x32f := (self, o, i, is) >>
     self.printf("vec_perm($1, $1, $2)", o.args),

));


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

NewRulesFor(InterpolateDFT, rec(
    InterpolateDFT_tSPL_PrunedDFT := rec(
       applicable := (self, nt) >> nt.hasTags() and ObjId(nt.params[2]) = fZeroPadMiddle,
       children := nt -> let(n:=nt.params[2].domain(), N := nt.params[1].range(), us := N/n, blk := n/2,
            [[
                TCompose([
                    Downsample(nt.params[1]),
                    PrunedDFT(N, 1, blk, [0, 2*us-1]),
                    TDiag(fConst(n, 1/n)),
                    DFT(n, -1)
                ]).withTags(nt.getTags())
            ]]),
       apply := (nt, C, cnt) -> C[1]
    )
));

_maxBkSize := function(exp)
    local i, l, d;

    if IsExp(exp) and ObjId(exp)=mul and IsValue(exp.args[1]) then
        return exp.args[1];
    fi;

    if IsInt(exp.eval()) or IsValue(exp.eval()) then return EvalScalar(exp); fi;
    l := Lambda(Filtered(exp.free(), IsLoopIndex), exp);
    d := List(spiral.sigma.GenerateData(l).tolist(), EvalScalar);
    return Gcd(d);
end;

_divideFunc := function(exp, dval)
    local l, d, i;

    if IsValue(exp) or IsInt(exp) then return exp/dval; fi;

    if IsExp(exp) and ObjId(exp)=mul and IsValue(exp.args[1]) then
        return ApplyFunc(mul, Concat([idiv(exp.args[1], dval)], Drop(exp.args, 1)));
    fi;

    if IsInt(exp.eval()) or IsValue(exp.eval()) then return EvalScalar(exp/dval); fi;
    i := Filtered(exp.free(), IsLoopIndex)[1];
    l := Lambda(i, exp);
    d := List(spiral.sigma.GenerateData(l).tolist(), EvalScalar)/dval;
    return FData(List(d, j->i.t.value(j))).at(i);
end;

NewRulesFor(Downsample, rec(
    Downsample_tag := rec(
        applicable := (self, nt) >> nt.hasTags(),
        apply := (nt, C, cnt) -> let(gcd := _maxBkSize(nt.params[1].domain()),
            NeedInterleavedComplex(ScatGath(fTensor(fId(_divideFunc(nt.params[1].domain(), gcd)), fId(gcd)), nt.params[1])))
    )
));

NewRulesFor(InterpolateSegmentDFT, rec(
    InterpolateSegmentDFT_tSPL_PrunedDFT := rec(
        applicable := (self, nt) >> nt.hasTags() and ObjId(nt.params[6]) = fZeroPadMiddle,
        children := nt -> let(tags := nt.getTags(),
                              numseg := nt.params[1], j := Ind(numseg), downsample := nt.params[5].at(j), upsample := nt.params[6],
                              insize := nt.params[3], overlap := nt.params[4], n := upsample.domain(), N := downsample.range(),
                              us := N/n, blk := n/2,
                              inp := n + (numseg-1)*(n-overlap), over := inp - insize, l := n - over, g := Gcd(n, l),

            [[
                Downsample(downsample).withTags(tags),
                PrunedDFT(N, 1, blk, [0, 2*us-1]).withTags(tags),
                DFT(n, -1).withTags(tags),
                PrunedDFT(n, -1, g, [0..l/g-1]).withTags(tags),
                InfoNt(j)
            ]]),
        apply := (nt, C, cnt) -> let(numseg := nt.params[1], _j := cnt[5].params[1], j := Ind(_j.range-1), downsample := nt.params[5].at(j), upsample := nt.params[6],
                              outsize := nt.params[2], insize := nt.params[3], overlap := nt.params[4], n := upsample.domain(), N := downsample.range(),
                              ids_rows := _sumSegDims(numseg-1, nt.params[5]),
#                              ids_rows := RulesStrengthReduce(outsize - nt.params[5].at(numseg-1).domain()),
                RowDirectSum(overlap,[
                    IRowDirSum(j, numseg-1, overlap,
                       SubstVars(Copy(C[1]), rec((_j.id) := j)) * C[2] * Diag(fConst(n, 1/n)) * C[3]).overrideDims([ids_rows, insize - C[4].dims()[2] + overlap]),
                    SubstVars(Copy(C[1]), rec((_j.id) := V(numseg-1))) * C[2] * Diag(fConst(n, 1/n)) * C[4]
                 ]).overrideDims([nt.params[2], nt.params[3]])
       )
    )
));


text:
# GAP has built-in support for quaternions

A := QuaternionAlgebra(Rationals);
# <algebra-with-one of dimension 4 over Rationals>

b := BasisVectors(Basis(A));
# [ e, i, j, k ]

q := [1, 2, 3, 4]*b;
# e+(2)*i+(3)*j+(4)*k

# Conjugate
ComplexConjugate(q);
# e+(-2)*i+(-3)*j+(-4)*k

# Division
1/q;
# (1/30)*e+(-1/15)*i+(-1/10)*j+(-2/15)*k

# Computing norm may be difficult, since the result would be in a quadratic field.
# Sqrt exists in GAP, but it is quite unusual: see ?E in GAP documentation, and the following example
Sqrt(5/3);
# 1/3*E(60)^7+1/3*E(60)^11-1/3*E(60)^19-1/3*E(60)^23-1/3*E(60)^31+1/3*E(60)^43-1/3*E(60)^47+1/3*E(60)^59

# However, the square of the norm is easy to compute
q*ComplexConjugate(q);
# (30)*e

q1 := [2, 3, 4, 5]*b;
# (2)*e+(3)*i+(4)*j+(5)*k

q2 := [3, 4, 5, 6]*b;
# (3)*e+(4)*i+(5)*j+(6)*k

q1*q2 - q2*q1;
# (-2)*i+(4)*j+(-2)*k

# Can't add directly to a rational, one must make a quaternion of it
r := 5/3*b[1];
# (5/3)*e
r + q;
# (8/3)*e+(2)*i+(3)*j+(4)*k

# For multiplication, no problem (we are in an algebra over rationals !)
r*q;
# (5/3)*e+(10/3)*i+(5)*j+(20/3)*k
5/3*q;
# (5/3)*e+(10/3)*i+(5)*j+(20/3)*k

# Negative
-q;
(-1)*e+(-2)*i+(-3)*j+(-4)*k

# While quaternions are built-in, you can define an algebra in GAP by specifying it's multiplication table.
# See tutorial, p. 60, and reference of the functions used below.

# A multiplication table of dimension 4.

T := EmptySCTable(4, 0);
SetEntrySCTable(T, 1, 1, [1, 1]);
SetEntrySCTable(T, 1, 2, [1, 2]);
SetEntrySCTable(T, 1, 3, [1, 3]);
SetEntrySCTable(T, 1, 4, [1, 4]);
SetEntrySCTable(T, 2, 1, [1, 2]);
SetEntrySCTable(T, 2, 2, [-1, 1]);
SetEntrySCTable(T, 2, 3, [1, 4]);
SetEntrySCTable(T, 2, 4, [-1, 3]);
SetEntrySCTable(T, 3, 1, [1, 3]);
SetEntrySCTable(T, 3, 2, [-1, 4]);
SetEntrySCTable(T, 3, 3, [-1, 1]);
SetEntrySCTable(T, 3, 4, [1, 2]);
SetEntrySCTable(T, 4, 1, [1, 4]);
SetEntrySCTable(T, 4, 2, [1, 3]);
SetEntrySCTable(T, 4, 3, [-1, 2]);
SetEntrySCTable(T, 4, 4, [-1, 1]);

A := AlgebraByStructureConstants(Rationals, T, ["e", "i", "j", "k"]);
b := GeneratorsOfAlgebra(A);

IsAssociative(A);
# true

IsCommutative(A);
# false

# Then, like above

q := [1, 2, 3, 4]*b;
# e+(2)*i+(3)*j+(4)*k

# However, as is, GAP does not know division or conjugate on this algebra.
# QuaternionAlgebra is useful as well for extensions of rationals,
# and this one _has_ conjugate and division, as seen previously.

# Try this on Q[z] where z is the square root of 5 (in GAP it's ER(5))
F := FieldByGenerators([ER(5)]);
A := QuaternionAlgebra(F);
b := GeneratorsOfAlgebra(A);

q := [1, 2, 3, 4]*b;
# e+(2)*i+(3)*j+(4)*k

# Conjugate and division

ComplexConjugate(q);
# e+(-2)*i+(-3)*j+(-4)*k

1/q;
# (1/30)*e+(-1/15)*i+(-1/10)*j+(-2/15)*k