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

# alpha is a function of P (divisor of N)
# m = N/P

# Frequency domain, u in [0, N)
# Per processor, total pts=2*m, total useful pts=m
# We start to see aliasing at 3*m/2, m=N/P.
W := (alpha, u) -> exp(-alpha * u^2);

# Time domain, t in [0, 1]
w := (alpha, t) -> sqrt(d_PI / alpha) * exp(-d_PI^2 * t^2 / alpha);

# digits - eg. 17, number of exact (unaffected by aliasing) decimal digits
gauss_alpha := (N, p, digits) -> let(m := N/p, 4/9 * digits * d_log(10) / m^2);

# how many samples of the filter to keep in the time domain, assuming gauss_alpha(N, p, digits)
# was used.
edge := digits -> IntDouble(d_ceil(2/(3*d_PI) * d_log(10) * digits));

Class(fPerExt, FuncClass, rec(
    def := (N, l, r) -> rec(N:=N, n:=N+l+r),
    lambda := self >> let(N := self.params[1], l:=self.params[2], r:=self.params[3],
 	j := Ind(self.n),
	Lambda(j, imod(j-l, N)))
));

# below we should not use HZ, since HZ assumes that n * stride = N, and if used as below
# assumption will be violated
#RewriteRules(RulesFuncSimp, rec(
#    fPerExt_H := ARule(fCompose, [@(1,fPerExt), @(2,H)], e ->
#	let(base := @(2).val.params[3], stride := @(2).val.params[4],
#	[ HZ(@(1).val.N, @(2).val.n, base - @(1).val.params[2], stride) ]))
#));

NewRulesFor(DFT, rec(

   DFT_Gauss := rec(
	forTransposition := true,
	minSize     := 4,
	maxSize     := false,
	p      := 4,
	digits := 15,

	applicable  := (self, t) >> (self.minSize<>false and Rows(t) > self.minSize) and
	                            not IsPrime(Rows(t)) and t.params[2] in [1,-1 mod Rows(t)],

        libApplicable := (self, t) >> eq(imod(t.params[1], self.p), 0),
	freedoms := (self, t) >> [],

	child := (self, t, freedoms) >> let(N := Rows(t), p := self.p,
	    [ DFT(2*N/p, t.params[2]), DFT(p, t.params[2]) ]),

	apply := meth(self, t,C,Nonterms)
	    local N, m, n, rot, k, p, alpha, time_filter, freq_filter, j;
	    N := Rows(t);
	    n := Rows(C[2]);
	    k := edge(self.digits);
	    p := self.p;
	    m := N/p;
	    alpha := gauss_alpha(N, p, self.digits);

	    time_filter := let(i:=Ind(2*k*p), Lambda(i, w(alpha, (i-k*p)/N)).setRange(TDouble));
	    freq_filter := let(i:=Ind(m),     Lambda(i, p/(2*N*W(alpha, i-m/2))).setRange(TDouble));

	    j := Ind(p);
	    return
	    Grp(
		IterDirectSum(j, j.range,
		    Diag(freq_filter) *
		    Gath(H(2*m, m, m*imod(j,2), 1)).toloop(4) * C[1]
		) *
		(Tr(2*N/p, p))
	    ) *
	    RowTensor(2*N/p, 2*k*p-p/2,
		C[2] *
		let(i:=Ind(p), Diag(Lambda(i, omega(2*p, i)).setRange(TComplex))) *
		Tensor(RowVec(List([0..2*k-1], x->(-1)^x)), I(p)) *
		Diag(time_filter)
	    ) *
	    Gath(fPerExt(N, k*p, k*p-p/2));
	end

   )
));

# scaling by p/(2*N) belongs to the time_filter, but freq_filter is shorter
#	    time_filter := List([-k*p .. k*p-1], i-> w(alpha, i/N));
#	    freq_filter := List([-m/2..m/2-1], i-> p/(2*N*W(alpha, i)));

#	    return
#	    Z(N, m/2) *
#	    Tensor(I(p/2),
#		   DirectSum(
#		       Diag(freq_filter) * Gath(fStack(H(2*m, m/2, 3*m/2, 1), H(2*m, m/2, 0, 1))) * C[1],
#		       Diag(freq_filter) * Gath(H(2*m, m, m/2, 1)) * C[1])
#	    ) *

testGauss := function(N, p, k)
    local s, me, them;
    DFT_Gauss.p := p;
    SwitchRulesByNameQuiet(DFT, [DFT_Gauss]);
    s := SPLRuleTree(ExpandSPL(DFT(N,k))[1]);
    me := MatSPL(s);
    them := MatSPL(DFT(N,k));
    VisualizeMat(me-them, " ");
    return [s, me, them];
end;

splGauss := function(N, p, k)
    local s, me, them;
    DFT_Gauss.p := p;
    SwitchRulesByNameQuiet(DFT, [DFT_Gauss]);
    s := SPLRuleTree(ExpandSPL(DFT(N,k))[1]);
    return s;
end;

sums := function(s, opts)
    local real;
    real := s.isReal();
    s := SumsSPL(s);
    s := SubstBottomUp(s, @.cond(IsNonTerminal), x->RecursStep(0,0,x));
    s := ApplyStrategy(s, opts.formulaStrategies.sigmaSpl, UntilDone);
    if not opts.generateComplexCode and not real then
        s := ApplyStrategy(RC(s), opts.formulaStrategies.rc, UntilDone); fi;
    return ApplyStrategy(s, opts.formulaStrategies.postProcess, UntilDone);
end;

# opts := CopyFields(SpiralDefaults, rec(generateComplexCode := true));
# opts := CopyFields(opts, rec(unparser := CMacroUnparser));

# N := 128;
# k := edge(17);
# p := 4;
# alpha := alpha(N, p, 17);

NewRulesFor(DFT, rec(

   DFT_PrunedGauss := rec(
	forTransposition := true,
	minSize     := 4,
	maxSize     := false,
	p := 2,
	digits := 15,

	applicable  := (self, t) >> (self.minSize<>false and Rows(t) > self.minSize) and
	                            not IsPrime(Rows(t)) and t.params[2] in [1,-1 mod Rows(t)],

	children := (self, t) >> let(N := Rows(t), p := self.p,
	    [[ DFT(2*N/p, t.params[2]), DFT3(p, t.params[2]) ]]), #DFT3

	apply := meth(self, t,C,Nonterms)
	    local N, m, n, mat, rot, k, p, alpha, time_filter, freq_filter, freq_filter2;
	    N := Rows(t);
	    n := Rows(C[2]);
            mat := MatSPL(C[2]);
	    k := edge(self.digits);
	    p := self.p;
	    m := N/p;
	    alpha := gauss_alpha(N, p, self.digits);
	    time_filter := List([-k*p .. k*p-1], i-> w(alpha, i/N).eval());
	    freq_filter := List([-m/2..m/2-1], i-> p/(2*N*W(alpha, i).eval()));

	    return
	    Diag(freq_filter) * Gath(H(2*m, m, 0, 1)) * C[1] *
	    RowTensor(2*N/p, 2*k*p-p/2,
		RowVec(mat[1]) *
		Tensor(RowVec(List([0..2*k-1], x->(-1)^x)), I(p)) *
		Diag(time_filter)) *
	    transforms.filtering.Ext(N, transforms.filtering.per(k * p),
		                        transforms.filtering.per(k * p - p/2));
	end
   )
));

#2N/p  * 2kp
#4Nk complex fmas overhead for filter (one stripe out), k=edge(digits)


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

mizuno:=[[1,2,3,4,5,6,7,
          9,10,11,12,
          13,14,16,17,
          19,18,20,21,
          23,24,25,27,
          26,28,30,31,
          32,33,35,34,
          37,38,39,40,
          41,44,45,46,
          48,49,51,52,
          53,55,57,58,
          59,61,63,64,
          66,69,70,71,
          75,76,80,82,
          85,89,93,97],
         [8,
          15,22,29,36,
          43,42,47,50,
          54,56,60,62,
          65,67,68,72,
          73,74,77,78,
          79,81,83,84,
          86,87,88,90,
          91,92,94,95,
          96,98,99,100,
          101,102,103,104,
          105,106,107,108,
          109,110,111,112,
          114,113,115,116,
          117,118,119,120]];

reps:=[[[1,1],[2,1],[3,1],[4,1],[5,1],[6,1],[7,1],[101,1]],          # E_8
       [[1,1],[2,1],[9,1],[10,1],[5,1],[6,1],[7,1],[101,1]],         # E_8(a_1)
       [[1,1],[3,1],[2,1],[9,1],[11,1],[12,1],[13,1],[101,1]],       # E_8(a_2)
       [[8,1],[9,1],[10,1],[11,1],[16,1],[12,1],[7,1],[101,1]],      # E_7+A_1
       [[1,1],[15,1],[16,1],[17,1],[6,1],[7,1],[101,1]],             # E_7
       [[8,1],[9,1],[10,1],[11,1],[16,1],[12,1],[13,1],[102,1]],     # D_8
       [[8,1],[15,1],[16,1],[17,1],[18,1],[19,1],[7,1],[101,1]],     # E_7(a_1)+A_1
       [[1,1],[15,1],[16,1],[17,1],[24,1],[13,1],[101,1]],           # E_7(a_1)
       [[1,1],[15,1],[16,1],[17,1],[18,1],[24,1],[13,1],[102,1]],    # D_8(a_1)       ?
       [[1,1],[15,1],[16,1],[17,1],[18,1],[19,1],[103,1]],           # D_7
       [[14,1],[15,1],[21,1],[17,1],[18,1],[29,1],[19,1],[101,1]],   # E_7(a_2)+A_1
       [[20,1],[21,1],[17,1],[23,1],[24,1],[7,1],[101,1]],           # E_7(a_2)
       [[14,1],[20,1],[21,1],[22,1],[18,1],[33,1],[13,1],[102,1]],   # A_8
       [[20,1],[21,1],[28,1],[23,1],[24,1],[7,1],[101,1]],           # E_6+A_1
       [[14,1],[22,1],[23,1],[24,1],[19,1],[25,1],[35,1],[102,1]],   # D_7(a_1)
#       [[14,1],[22,1],[23,1],[24,1],[19,1],[25,1],[102,1]],          # (D_7(a_1))_2
       [[20,1],[21,1],[23,1],[24,1],[7,1],[101,1]],                  # E_6
       [[20,1],[21,1],[22,1],[23,1],[24,1],[25,1],[102,1],[103,1]],  # D_8(a_2)
       [[8,1],[28,1],[23,1],[24,1],[25,1],[30,1],[102,1]],           # D_6+A_1
       [[20,1],[21,1],[28,1],[23,1],[24,1],[10,1],[40,1],[103,1]],   # A_7
       [[20,1],[21,1],[28,1],[23,1],[24,1],[19,1],[103,1]],          # (A_7)_3
       [[20,1],[21,1],[28,1],[23,1],[31,1],[101,1],[102,1]],         # E_6(a_1)+A_1
       [[1,1],[28,1],[29,1],[30,1],[31,1],[101,1]],                  # D_6
       [[20,1],[21,1],[28,1],[29,1],[25,1],[31,1],[103,1]],          # D_7(a_2)
       [[20,1],[21,1],[23,1],[31,1],[101,1],[102,1]],                # E_6(a_1)
       [[26,1],[28,1],[27,1],[29,1],[30,1],[31,1],[102,1],[103,1]],  # D_5+A_2
#       [[26,1],[28,1],[27,1],[29,1],[30,1],[31,1],[102,1]],          # (D_5+A_2)_2
       [[20,1],[21,1],[28,1],[29,1],[30,1],[31,1],[105,1]],          # A_6+A_1
       [[28,1],[27,1],[29,1],[30,1],[31,1],[39,1],[102,1]],          # D_6(a_1)+A_1
       [[28,1],[29,1],[33,1],[31,1],[34,1],[101,1]],                 # D_6(a_1)
       [[29,1],[21,1],[20,1],[30,1],[31,1],[105,1]],                 # A_6
#       [[32,1],[33,1],[30,1],[32,1],[40,1],[101,1]],                 # D_5+A_1        ??? two times 32
       [[26,1],[29,1],[35,1],[38,1],[30,1],[44,1],[47,1],[105,1]],   # 2A_4           ? 26
       [[32,1],[37,1],[29,1],[31,1],[34,1],[105,1],[107,1]],         # A_5+A_2
       [[32,1],[37,1],[27,1],[33,1],[40,1],[36,1],[105,1]],          # A_5+2A_1
       [[32,1],[29,1],[31,1],[34,1],[106,1],[108,1]],                # D_6(a_2)
       [[30,1],[31,1],[32,1],[33,1],[101,1]],                        # D_5
       [[37,1],[33,1],[40,1],[34,1],[36,1],[41,1],[105,1]],          # D_5(a_1)+A_2   ? 26/20
       [[26,1],[27,1],[40,1],[41,1],[106,1],[107,1]],                # (A_5+A_1)'
       [[32,1],[27,1],[37,1],[33,1],[36,1],[105,1]],                 # (A_5+A_1)''
       [[37,1],[38,1],[40,1],[49,1],[41,1],[104,1],[108,1]],         # D_4+A_2
 #      [[37,1],[38,1],[40,1],[49,1],[41,1],[104,1]],                 # (D_4+A_2)_2
       [[37,1],[38,1],[40,1],[34,1],[36,1],[106,1],[107,1]],         # A_4+A_3
       [[20,1],[21,1],[49,1],[106,1],[107,1]],                       # A_5
       [[37,1],[42,1],[43,1],[44,1],[45,1],[103,1]],                 # D_5(a_1)+A_1
       [[37,1],[38,1],[39,1],[40,1],[41,1],[108,1],[109,1]],         # A_4+A_2+A_1

       [[38,1],[40,1],[39,1],[41,1],[108,1],[109,1]],       # A_4+A_2
       [[42,1],[43,1],[48,1],[49,1],[102,1]],               # D_5(a_1)
       [[47,1],[48,1],[49,1],[53,1],[101,1]],               # D_4+A_1
       [[42,1],[43,1],[44,1],[45,1],[107,1],[108,1]],       # A_4+2A_1
       [[47,1],[48,1],[49,1],[101,1]],                      # D_4
       [[42,1],[48,1],[49,1],[106,1],[107,1]],              # A_4+A_1
       [[37,1],[38,1],[39,1],[49,1],[108,1],[113,1]],       # 2A_3
       [[42,1],[48,1],[106,1],[107,1]],                     # A_4
       [[42,1],[44,1],[48,1],[49,1],[110,1],[113,1]],       # D_4(a_1)+A_2
       [[46,1],[47,1],[48,1],[49,1],[112,1],[113,1]],       # A_3+A_2+A_1
       [[46,1],[47,1],[48,1],[112,1],[116,1],[119,1]],      # A_3+A_2
#       [[46,1],[47,1],[48,1],[112,1],[116,1]],              # (A_3+A_2)_2
       [[53,1],[47,1],[48,1],[115,1],[116,1]],              # D_4(a_1)+A_1
       [[53,1],[54,1],[55,1],[112,1],[113,1]],              # A_3+2A_1
       [[47,1],[48,1],[115,1],[116,1]],                     # D_4(a_1)
       [[53,1],[54,1],[55,1],[117,1],[118,1],[119,1]],      # 2A_2+2A_1
       [[55,1],[56,1],[112,1],[113,1]],                     # A_3+A_1
       [[56,1],[57,1],[117,1],[118,1],[119,1]],             # 2A_2+A_1
       [[63,1],[106,1],[107,1]],                            # A_3
       [[56,1],[57,1],[117,1],[118,1]],                     # 2A_2
       [[60,1],[126,1],[127,1],[128,1]],                    # A_2+2A_1
       [[58,1],[59,1],[123,1],[124,1],[125,1]],             # A_2+3A_1
       [[63,1],[128,1],[129,1]],                            # A_2+A_1
       [[127,1],[130,1]],                                   # A_2
       [[63,1],[135,1],[136,1],[137,1]],                    # 4A_1
       [[141,1],[142,1],[143,1]],                           # 3A_1
       [[150,1],[151,1]],                                   # 2A_1
       [[157,1]]                                            # A_1
       ];

labels:=["E_8",
         "E_8(a_1)",
         "E_8(a_2)",
         "E_7+A_1",
         "E_7",
         "D_8",
         "E_7(a_1)+A_1",
         "E_7(a_1)",
         "D_8(a_1)",
         "D_7",
         "E_7(a_2)+A_1",
         "E_7(a_2)",
         "A_8",
         "E_6+A_1",
         "D_7(a_1)",
#         "(D_7(a_1))_2",
         "E_6",
         "D_8(a_2)",
         "D_6+A_1",
         "A_7",
         "(A_7)_3",
         "E_6(a_1)+A_1",
         "D_6",
         "D_7(a_2)",
         "E_6(a_1)",
         "D_5+A_2",
#         "(D_5+A_2)_2",
         "A_6+A_1",
         "D_6(a_1)+A_1",
         "D_6(a_1)",
         "A_6",
#         "D_5+A_1",
         "2A_4",
         "A_5+A_2",
         "A_5+2A_1",
         "D_6(a_2)",
         "D_5",
         "D_5(a_1)+A_2",
         "(A_5+A_1)'",
         "(A_5+A_1)''",
         "D_4+A_2",
#         "(D_4+A_2)_2",
         "A_4+A_3",
         "A_5",
         "D_5(a_1)+A_1",
         "A_4+A_2+A_1",

         "A_4+A_2",
         "D_5(a_1)",
         "D_4+A_1",
         "A_4+2A_1",
         "D_4",
         "A_4+A_1",
         "2A_3",
         "A_4",
         "D_4(a_1)+A_2",
         "A_3+A_2+A_1",
         "A_3+A_2",
#         "(A_3+A_2)_2",
         "D_4(a_1)+A_1",
         "A_3+2A_1",
         "D_4(a_1)",
         "2A_2+2A_1",
         "A_3+A_1",
         "2A_2+A_1",
         "A_3",
         "2A_2",
         "A_2+2A_1",
         "A_2+3A_1",
         "A_2+A_1",
         "A_2",
         "4A_1",
         "3A_1",
         "2A_1",
         "A_1",
         ];

I:=[ [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ] ];#, [  ], [  ], [  ], [  ], [  ] ];

II:=[ [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ] ];#, [  ], [  ], [  ], [  ], [  ] ];

tori:=[ [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ], [  ] ];#, [  ], [  ], [  ], [  ], [  ] ];

#blocks:=[];
blocks:=[[[55, 1], [45, 1], [39, 1], [27, 3], [19, 1], [9, 1]], [[27, 8], [25, 1], [7, 1]], [[27, 6], [25, 1], [21, 1], [19, 1], [9, 2], [3, 1]], [[27, 4], [25, 1], [19, 1], [18, 2], [15, 1], [12, 2], [9, 2], [3, 1]], [[27, 3], [26, 2], [19, 1], [18, 2], [15, 1], [12, 2], [9, 2], [1, 3]], [[27, 2], [25, 1], [21, 2], [18, 2], [15, 2], [12, 2], [9, 3], [7, 1], [3, 1]], [[27, 1], [21, 1], [19, 2], [18, 4], [15, 1], [9, 8], [3, 1]], [[27, 1], [20, 2], [19, 1], [18, 4], [15, 1], [9, 8], [1, 3]], [[19, 1], [15, 8], [9, 9], [3, 8], [2, 2]], [[23, 1], [22, 2], [19, 1], [16, 2], [15, 1], [13, 3], [12, 2], [11, 1], [10, 2], [7, 1], [4, 2], [3, 1], [1, 3]], [[19, 1], [18, 4], [15, 1], [12, 2], [9, 12], [7, 1], [3, 1]], [[19, 1], [18, 4], [15, 1], [12, 2], [9, 11], [8, 2], [1, 3]], [[19, 1], [18, 2], [15, 3], [13, 3], [9, 9], [7, 1], [6, 2], [3, 3]], [[19, 1], [15, 8], [9, 9], [3, 8], [2, 2]], [[19, 1], [17, 2], [15, 3], [12, 2], [11, 4], [9, 5], [6, 4], [3, 4], [1, 1]], [[19, 1], [15, 8], [9, 9], [3, 7], [1, 7]], [[9, 25], [8, 2], [2, 2], [1, 3]], [[19, 1], [17, 1], [16, 2], [15, 2], [12, 2], [11, 2], [10, 4], [9, 3], [6, 4], [3, 2], [2, 2], [1, 3]], [[9, 3], [7, 5], [6, 16], [5, 3], [3, 24], [1, 3]], [[9, 25], [8, 2], [2, 2], [1, 3]], [[9, 26], [7, 1], [3, 1], [2, 2]], [[19, 1], [16, 4], [15, 1], [11, 5], [10, 4], [9, 2], [6, 4], [3, 1], [1, 10]], [[9, 23], [7, 1], [6, 2], [5, 2], [4, 2], [3, 1], [1, 1]], [[9, 26], [7, 1], [1, 7]], [[9, 21], [7, 2], [6, 4], [3, 7]], [[9, 19], [8, 2], [7, 4], [6, 2], [3, 6], [1, 3]], [[9, 20], [8, 2], [7, 1], [6, 4], [3, 6], [1, 3]], [[9, 19], [8, 4], [6, 4], [3, 5], [1, 6]], [[9, 19], [7, 8], [3, 5], [1, 6]], [[9, 14], [7, 2], [6, 8], [5, 6], [3, 10]], [[9, 13], [8, 2], [7, 1], [6, 8], [5, 3], [4, 6], [3, 6], [1, 3]], [[9, 10], [7, 1], [6, 16], [3, 17], [2, 2]], [[9, 13], [7, 5], [6, 4], [5, 5], [4, 8], [3, 3], [1, 6]], [[9, 19], [7, 8], [1, 21]], [[9, 9], [8, 2], [6, 16], [3, 17], [2, 2]], [[9, 9], [8, 2], [6, 16], [3, 16], [2, 2], [1, 3]], [[9, 10], [7, 1], [6, 16], [3, 16], [1, 7]], [[9, 3], [8, 8], [7, 6], [6, 10], [3, 8], [2, 8], [1, 15]], [[9, 7], [8, 2], [7, 6], [6, 6], [5, 6], [4, 6], [3, 10], [2, 2], [1, 3]], [[9, 9], [8, 2], [6, 16], [3, 15], [1, 10]], [[9, 5], [8, 6], [7, 4], [6, 12], [3, 11], [2, 8], [1, 6]], [[9, 3], [8, 2], [7, 1], [6, 18], [5, 3], [3, 25]], [[9, 3], [7, 5], [6, 16], [5, 3], [3, 24], [1, 3]], [[9, 3], [8, 8], [7, 6], [6, 10], [3, 8], [2, 8], [1, 15]], [[9, 2], [8, 6], [7, 13], [6, 6], [3, 2], [2, 14], [1, 21]], [[9, 1], [8, 4], [7, 3], [6, 12], [5, 7], [4, 4], [3, 17], [2, 4], [1, 4]], [[9, 2], [7, 25], [3, 1], [1, 52]], [[9, 1], [8, 2], [7, 7], [6, 8], [5, 9], [4, 6], [3, 12], [2, 6], [1, 9]], [[8, 4], [7, 6], [6, 4], [5, 10], [4, 16], [3, 6], [2, 4], [1, 10]], [[9, 1], [7, 11], [5, 21], [3, 11], [1, 24]], [[7, 1], [6, 14], [5, 6], [3, 42], [1, 1]], [[7, 1], [6, 14], [5, 3], [4, 6], [3, 37], [2, 2], [1, 3]], [[7, 1], [6, 12], [5, 5], [4, 8], [3, 34], [1, 10]], [[7, 1], [6, 8], [5, 12], [4, 6], [3, 28], [2, 8], [1, 9]], [[7, 1], [6, 6], [5, 11], [4, 14], [3, 19], [2, 12], [1, 13]], [[7, 1], [6, 2], [5, 24], [3, 27], [1, 28]], [[3, 80], [2, 4]], [[7, 1], [6, 2], [5, 15], [4, 18], [3, 10], [2, 14], [1, 24]], [[3, 79], [2, 2], [1, 7]], [[7, 1], [5, 11], [4, 32], [3, 1], [1, 55]], [[3, 78], [1, 14]], [[3, 65], [2, 16], [1, 21]], [[3, 70], [2, 14], [1, 10]], [[3, 58], [2, 20], [1, 34]], [[3, 57], [1, 77]], [[3, 44], [2, 40], [1, 36]], [[3, 31], [2, 50], [1, 55]], [[3, 14], [2, 64], [1, 78]], [[3, 1], [2, 56], [1, 133]]];

#ords:=[];
ords:=[ 81, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 27, 9, 27, 9, 9, 9, 27, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 3, 9, 3, 9, 3, 3, 3, 3, 3, 3, 3, 3, 3 ];

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

ext: gi


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

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

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

highest_weight:=[0,2];# corresponding highest weight module is the Lie algebra

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

#
# p>2
#

getRep:=function(k,t)
    local u,ua,ub,uab,uabb;
    #ua:=ext_root_group(k,1,t);
    ub:=ext_root_group(k,2,t);
    uab:=ext_root_group(k,3,t);
    uabb:=ext_root_group(k,4,t^2);
    u:=ub*uab*uabb;

    return u;
end;

handle:=function(k)
    local id_mat,u,e,h,h1,h2,ese,evh,esh,esh0,index_ev0,si,inter;
    u:=getRep(k,1);
    id_mat:=u^0;
    e:=u-id_mat;
    ese:=Eigenspaces(Rationals,e)[1]; #eigenvalues 0, centralizer

    # torus
    h1:=ext_e(k,cb0[7]);
    h2:=ext_e(k,cb0[8]);;
    h:=h1+h2;
    esh:=Eigenspaces(Rationals,h);
    evh:=Eigenvalues(Rationals,h);

    # centralizer of torus
    esh0:=[];
    index_ev0:=Position(evh,0);
    esh0:=esh[index_ev0];

    #
    # sum and intersection
    #
    si:=SumIntersectionMat(Basis(ese),Basis(esh0));
    inter:=si[2];
    return inter;
end;

max_stars_repo_name: iuliansimion/groupsSB

max_issues_repo_name: iuliansimion/groupsSB

max_forks_repo_name: iuliansimion/groupsSB


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

#
# NB: Contents of this file should be moved to a directory local to the
#     particular Spiral Extension
#

Class(f2DTr, Tr, rec(
    advdomain := self >> [[self.params[2], self.params[1]]],
    advrange := self >> [[self.params[1], self.params[2]]]
));

Class(f2DTrExplicitL, f2DTr);
Class(f2DTrExplicitR, f2DTr);

Class(f4DHalfFlip, PermClass, rec(
    def := (a1,a2,b1,b2) -> rec(size := a1*a2*b1*b2),
    lambda := self >> let(
        a1 := self.params[1], a2 := self.params[2],
        b1 := self.params[3], b2 := self.params[4],
        fTensor(Tr(a1, a2), Tr(b1, b2)).lambda()
        ),
    domain := self >> self.params[1] * self.params[2] * self.params[3] * self.params[4],
    range := self >> self.params[1] * self.params[2] * self.params[3] * self.params[4],

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

    transpose := self >> Error("not supported"),

    advdomain := self >> [[self.params[4], self.params[3], self.params[2], self.params[1]]],
    advrange := self >> [[self.params[3], self.params[4], self.params[1], self.params[2]]]
));

Class(f4DMacroFlip, PermClass, rec(
    def := (a1,a2,b1,b2) -> rec(size := a1*a2*b1*b2),
    lambda := self >> let(
        a1 := self.params[1], a2 := self.params[2],
        b1 := self.params[3], b2 := self.params[4],
        Tr(a1*a2, b1*b2).lambda()),
#         i := Ind(a1*a2*b1*b2),
#         k := a1*a2, str := b1*b2,
#         Lambda(i, add(mul(b2, idiv(idiv(i, k), b2)), imod(idiv(i, k), b2))
#             + str * add(mul(a2, idiv(imod(i, k), a2)), imod(imod(i, k), a2)))),

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

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

    transpose := self >> Error("not supported"),

    advdomain := self >> [[self.params[3], self.params[4], self.params[1], self.params[2]]],
    advrange := self >> [[self.params[1], self.params[2], self.params[3], self.params[4]]]
));

Class(f4DFullFlip, PermClass, rec(
    def := (a1,a2,b1,b2) -> rec(size := a1*a2*b1*b2),
    lambda := self >> let(
        a1 := self.params[1], a2 := self.params[2],
        b1 := self.params[3], b2 := self.params[4],
        i := Ind(a1*a2*b1*b2),
        Lambda(i, imod(i,b2) * a1 * a2 * b1
            + imod(idiv(i,b2),b1)*a1*a2
            + imod(idiv(i, b2*b1),a2) * a1
            + imod(idiv(i, b2*b1*a2), a1))),
    domain := self >> self.params[1] * self.params[2] * self.params[3] * self.params[4],
    range := self >> self.params[1] * self.params[2] * self.params[3] * self.params[4],

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

    transpose := self >> Error("not supported"),

    advdomain := self >> [[self.params[4], self.params[3], self.params[2], self.params[1]]],
    advrange := self >> [[self.params[1], self.params[2], self.params[3], self.params[4]]]
));

#F
#F Multi-dimensional fTensor's for OL
#F

Class(f2DTensor, fTensor, rec(
    advdomain := self >> [List(self._children, x->x.domain())],
    advrange := self >> [List(self._children, x->x.range())]
));

Class(f4DTensor, fTensor, rec(
    advdomain := self >> [List(self._children, x->x.domain())],
    advrange := self >> [List(self._children, x->x.range())]
));


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

avarnames:=List([1..7],i->Concatenation("a_{",String(i),"}"));
PR:=PolynomialRing(GF(5),avarnames);

Print("Created PolynomialRing ",PR,"\n");
Print("\t[",IsPolynomialRing(PR)=true,"] IsPolynomialRing(PR)=true\n");


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

NewRulesFor(DFT, rec(
	DFT_tSPL_CT_Mem := rec(
		switch := true,

		applicable := (self, nt) >> nt.params[1] > 2
			and not IsPrime(nt.params[1])
			and nt.hasTag(AMem),

		children := nt -> let(
			memtag := nt.getTag(AMem),
			Map2(
				Filtered(DivisorPairs(nt.params[1]), (m,n) -> n mod memtag.block = 0),
				(m,n) -> [
					TCompose([
						TGrp(TCompose([
							TTensorI(DFT(m, nt.params[2] mod n), n, AVec, AVec),
							TDiag(fPrecompute(Tw1(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]
	)
));


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

Class(spu_common, SIMD_VMX, rec(
    file := "cell",
    unparser := CellUnparser,
    #compileStrategy := self >> IndicesCS2_FMA,
    compileStrategy := self >> IndicesCS2,
    useDeref := true,

    vadd    := "spu_add",
    vsub    := "spu_sub",
    vmul    := "spu_mul",
    vlshift := "spu_sl",
    vrshift := "spu_sr",
    vbinand := "spu_and",

    autolib := rec(
       #includes := () -> ["<vec_types.h>", "<altivec.h>", "<spu2vmx.h>" ], # This is for PowerPC (PPE)
       includes := () -> ["<spu_intrinsics.h>", "<include/sp_general_malloc.h>",
                   "<simdmath.h>", "<simdmath/negated2.h>", "<simdmath/negatef4.h>"],
       timerIncludes := () -> ["<include/sp_spe_timer.h>" ],
    ),

    backendConfig := rec(
                        profile := default_profiles.linux_cellmultiSPU_gcc,
                        measureFunction := _StandardMeasureVerify
                     ),
));

Class(spu_4x32f, spu_common, rec(
    active := true,
    isFixedPoint := false,
    info := "SPU 4 x 32-bit float",
    v := 4,
    t := TVect(TReal, 4),
    globalUnrolling := 256,
    ctype := "float",
    vtype := "vector float",
    stype := "__attribute__ ((aligned(16))) float",
    instr := [ vperm_4x32f_spu, vuperm_4x32f_spu ],
    bits := 32,
    isFloat := true,
    isFix := false,
    includes := () -> ["<spu_intrinsics.h>", "<omega32.h>", "<simdmath.h>", "<simdmath/negated2.h>", "<simdmath/negatef4.h>"],
    infix_op := false,
    infix_assign := true,
    vmuladd := false,
    dupload := (y, x) -> assign(y, vuperm_4x32f_spu(promote_spu4x32f(x, 0), [1,1,1,1])),
    #duploadn takes a vector and replicates scalar number n in a vector
    duploadn := (y, x, n) -> assign(y, vuperm_4x32f_spu(x, [n,n,n,n])),
    vzero   := vzero_4x32f,
    vconst := vconstpr_av,
    vconstv := "(vector float)",
    vconst1 := "(vector float)spu_splats",
    splopts := rec(precision := "single"),
    countrec := rec( ops := [[add, sub], [mul], [fma, fms, nfma], [vperm_4x32f_spu]],
                     printstrings := ["[adds]", "[mults]", "[fmas]", "[vperms]"],
                     type := "TVect",
                     arithcost := (self, opcount) >> opcount[1]+opcount[2]+(2*opcount[3])
                     ),

    #NOTE:
    # - If operating on a NAN has no penalty, we shouldn't bother promoting and
    # shuffling to zero, for instance, where we could simply load directly.

    #NOTE:
    # - Find out if these are the fastest/most efficient (lowest register
    #   usage). For instance, how does a
    #   (promote,promote,perm,promote,promote,perm,perm) compare to a (promote,promote,perm,insert,insert)?
    # - Also, should we use temp variables here or not? (might help compiler propagation passes?)
    # - How to reuse loadc code here?

    #----------------------------------------
    #NOTE: Find out if ObjId(x)=nth should hold for all conditions. Doesn't seem to.
    #   (y,x) -> Checked(ObjId(x)=nth, assign(y, vperm_4x32f_spu(promote_spu4x32f(x, 0), vzero_spu(), [1, 128, 128, 128]))),
    svload := [[ # load using subvectors of length 1
              (y,x,opts) -> assign(y, vuperm_4x32f_spu(promote_spu4x32f(x[1], 0), [1, 128, 128, 128])),
              (y,x,opts) -> assign(y, vperm_4x32f_spu(promote_spu4x32f(x[1], 0), promote_spu4x32f(x[2], 0), [1,   5, 128, 128])),
              (y,x,opts) -> assign(y, insert_spu4x32f(x[3], vperm_4x32f_spu(promote_spu4x32f(x[1], 0), promote_spu4x32f(x[2], 0), [1,   5, 128, 128]), 2)),
              (y,x,opts) -> assign(y, insert_spu4x32f(x[4], insert_spu4x32f(x[3], insert_spu4x32f(x[2], promote_spu4x32f(x[1], 0), 1), 2), 3) ),
              ],
              [ # load using subvectors of length 2
              (y,x,opts) -> assign(y, vperm_4x32f_spu(promote_spu4x32f(x[1], 0), promote_spu4x32f(nth(x[1].loc, x[1].idx+1), 0), [1, 5, 128, 128])),
              (y,x,opts) -> assign(y, insert_spu4x32f(nth(x[2].loc, x[2].idx+1), insert_spu4x32f(x[2], insert_spu4x32f(nth(x[1].loc, x[1].idx+1), promote_spu4x32f(x[1], 0), 1), 2), 3) ),
              ]],
    #----------------------------------------
    #NOTE: loadu3 and loadu4 need work, prolly mostly in the unparser.
#     loadc := [      # load contiguous unaligned
#        (y,x) -> Checked(ObjId(x)=nth, assign(y, vuperm_4x32f_spu(promote_spu4x32f(x, 0), [1, 128, 128, 128]))),
#        (y,x) -> Checked(ObjId(x)=nth, assign(y, vperm_4x32f_spu(promote_spu4x32f(x, 0), promote_spu4x32f(nth(x.loc, x.idx+1), 0), [1, 5, 128, 128]))),

# #NOTE: This might be the best way to do it. Need to get it to work (past the ptr, ptr+1 problems).
# #(y,x) -> assign(y, vperm_4x32f_spu(bin_or(slqwbyte_spu4x32f(x.toPtr(TVect(TDouble, 4))), rlmaskqwbyte_spu4x32f(x.toPtr(TVect(TDouble, 4)))), vzero_spu(), [1,2,3,128])),
# #(y,x) -> assign(y, bin_or(slqwbyte_spu4x32f(x.toPtr(TVect(TDouble, 4))), rlmaskqwbyte_spu4x32f(x.toPtr(TVect(TDouble, 4))))),

# #NOTE: This is kind of a go-between for now. This bad because it's hacked via the unparser, but good because it only does 2 mem references
#        (y,x) -> Checked(ObjId(x)=nth, assign(y, vuperm_4x32f_spu(vloadu4_spu4x32f(x), [1, 2, 3, 128]))),
#        (y,x) -> Checked(ObjId(x)=nth, assign(y, vloadu4_spu4x32f(x)))
#     ],

###########################################################################
###### Vas, please check
###########################################################################

    # keep the n lower scalars and zero the other ones
    optional_mask :=  (c, n, opts) -> When(IsBound(opts.trueSVSemantics) and opts.trueSVSemantics and not(n=4),
        vuperm_4x32f_spu(c, List([1..4],x->When(x<=n, x, 128))),
        c),

    # load contiguous with unaligned loads
    loadc := (self, sv, opts) >> Cond(sv=1,
        ((y,x) -> Checked(ObjId(x)=nth, assign(y, vuperm_4x32f_spu(promote_spu4x32f(x, 0), [1, 128, 128, 128])))),
        sv=2,
        (y,x) -> Checked(ObjId(x)=nth, assign(y, vperm_4x32f_spu(promote_spu4x32f(x, 0), promote_spu4x32f(nth(x.loc, x.idx+1), 0), [1, 5, 128, 128]))),
#NOTE: loadu3 and loadu4 need work, prolly mostly in the unparser.
#NOTE: This is kind of a go-between for now. This bad because it's hacked via the unparser, but good because it only does 2 mem references
        sv=3,
        (y,x) -> Checked(ObjId(x)=nth, assign(y, vuperm_4x32f_spu(vloadu4_spu4x32f(x), [1, 2, 3, 128]))),
        sv=4,
        (y,x) -> Checked(ObjId(x)=nth, assign(y, vloadu4_spu4x32f(x))),
        Error("bad arguments for loadc")),

    #load contiguous + known alignment -> using 2 aligned load to be smarter
    loadc_align := (self, sv, align, opts)>>
    ((y,x,addr) -> let(
        v1 := nth(nth(x,add(addr,-align)).toPtr(TVect(TReal, 4)),0),
        v2 := nth(nth(x,add(addr,4-align)).toPtr(TVect(TReal, 4)), 0),
        m := x-> self.optional_mask(x, sv, opts),
        Cond(align=0,
            assign(y, m(v1)),
            assign(y, vperm_4x32f_spu(v1, v2, List([1..4], x-> When(x<=sv, align+x, 128))))
            ))),

###########################################################################
###### </Vas>
###########################################################################

    #----------------------------------------
    svstore := [[ # store using subvectors of length 1
            (y,x,opts) ->        assign(y[1], extract_spu4x32f(x, 0)),

            (y,x,opts) -> chain( assign(y[1], extract_spu4x32f(x, 0)),
                            assign(y[2], extract_spu4x32f(x, 1)) ),

            (y,x,opts) -> chain( assign(y[1], extract_spu4x32f(x, 0)),
                            assign(y[2], extract_spu4x32f(x, 1)),
                            assign(y[3], extract_spu4x32f(x, 2)) ),

            (y,x,opts) -> chain( assign(y[1], extract_spu4x32f(x, 0)),
                            assign(y[2], extract_spu4x32f(x, 1)),
                            assign(y[3], extract_spu4x32f(x, 2)),
                            assign(y[4], extract_spu4x32f(x, 3)) ),
              ],

              [ # store using subvectors of length 2
            (y,x,opts) -> chain( assign(y[1],                      extract_spu4x32f(x, 0)),
                            assign(nth(y[1].loc, y[1].idx+1), extract_spu4x32f(x, 1)) ),

            (y,x,opts) -> chain( assign(y[1],                      extract_spu4x32f(x, 0)),
                            assign(nth(y[1].loc, y[1].idx+1), extract_spu4x32f(x, 1)),
                            assign(y[2],                      extract_spu4x32f(x, 2)),
                            assign(nth(y[2].loc, y[2].idx+1), extract_spu4x32f(x, 3)) ),
              ]],
    #----------------------------------------
    storec := [    # store contiguous unaligned
            (y,x) ->        assign(y, extract_spu4x32f(x, 0)),
            (y,x) -> chain( assign(y,                   extract_spu4x32f(x, 0)),
                            assign(nth(y.loc, y.idx+1), extract_spu4x32f(x, 1)) ),
            (y,x) -> chain( assign(y,                   extract_spu4x32f(x, 0)),
                            assign(nth(y.loc, y.idx+1), extract_spu4x32f(x, 1)),
                            assign(nth(y.loc, y.idx+2), extract_spu4x32f(x, 2)) ),
            (y,x) -> chain( assign(y,                   extract_spu4x32f(x, 0)),
                            assign(nth(y.loc, y.idx+1), extract_spu4x32f(x, 1)),
                            assign(nth(y.loc, y.idx+2), extract_spu4x32f(x, 2)),
                            assign(nth(y.loc, y.idx+3), extract_spu4x32f(x, 3)) ),
        ],
    reverse := (y,x) -> assign(vref(y,0,4), vuperm_4x32f_spu(vref(x,0,4), [4, 3, 2, 1])),

    #----------------------------------------
    # support for VS and VS.transpose()

    shl1 :=  (y,x,opts) -> assign(y, vuperm_4x32f_spu(x, [128,1,2,3])),
    shl2 :=  (y,x,opts) -> assign(y, vperm_4x32f_spu(x[1], x[2], [4,5,6,7])),
    shr1 :=  (y,x,opts) -> assign(y, vuperm_4x32f_spu(x, [2, 3, 4, 128])),
    shr2 :=  (y,x,opts) -> assign(y, vperm_4x32f_spu(x[1], x[2], [2,3,4,5])),

    # support for VO1dsJ(n, v)
    shrev := (y,x,opts) -> assign(y, vperm_4x32f_spu(x[1], x[2], [1, 8, 7, 6])),

    swap_cx := (y, x, opts) -> assign(y, vuperm_2x64f_spu(x, [2,1,4,3])),
    RCVIxJ2 := (y, x, opts) -> assign(y, vuperm_2x64f_spu(x, [2,1,4,3]))

));

Class(spu_2x64f, spu_common, rec(
    active  := true,
    isFixedPoint := false,
    info    := "SPU 2 x 64-bit double",
    v       := 2,
    t := TVect(TReal, 2),
    globalUnrolling := 256,
    ctype   := "double",
    vtype   := "vector double",
    stype   := "__attribute__ ((aligned(16))) double",
    instr   := [ vperm_2x64f_spu, vuperm_2x64f_spu ],
    bits    := 64,
    isFloat := true,
    isFix := false,
    includes := () -> ["<spu_intrinsics.h>", "<omega32.h>"],
    infix_op := false,
    infix_assign := true,
    vmuladd := false,

    #NOTE: fix these
    #dupload := (y, x) -> assign(y, vuperm_4x32f_spu(promote_spu4x32f(x, 0), [1,1,1,1])),
    #duploadn takes a vector and replicates scalar number n in a vector
    #uploadn := (y, x, n) -> assign(y, vuperm_4x32f_spu(x, [n,n,n,n])),
    vzero   := vzero_2x64f,

    vconstv := "(vector double)",
    vconst1 := "(vector double)spu_splats",
    splopts := rec(precision := "double"),

    #NOTE: fix these
    #upload := (y, x) -> assign(y, vperm_2x64f_spu(promote_spu2x64f(x, 0), vzero_2x64f(), [1,1])),
    #uploadn := (y, x, n) -> assign(y, vuperm_2x64f_spu(x, [n,n])),

    # NOTE: Check if each of these are correct/efficient.
    svload := [[ # load using subvectors of length 1
              (y,x,opts) -> assign(y, vuperm_2x64f_spu(promote_spu2x64f(x[1], 0), [1, 128])),
              (y,x,opts) -> assign(y, vperm_2x64f_spu(promote_spu2x64f(x[1], 0), promote_spu2x64f(x[2], 0), [1,   3])),
              ],
              [ # load using subvectors of length 2
              (y,x,opts) -> assign(y, vperm_2x64f_spu(promote_spu2x64f(x[1], 0), promote_spu2x64f(nth(x[1].loc, x[1].idx+1), 0), [1, 3])),
              ]],

#     loadc := [      # load contiguous unaligned
#        (y,x) -> Checked(ObjId(x)=nth, assign(y, vuperm_2x64f_spu(promote_spu2x64f(x, 0), [1, 128]))),
#        (y,x) -> Checked(ObjId(x)=nth, assign(y, vperm_2x64f_spu(promote_spu2x64f(x, 0), promote_spu2x64f(nth(x.loc, x.idx+1), 0), [1, 3]))),
#        ],

###########################################################################
###### NOTE: Vas, please check
###########################################################################

    # keep the n lower scalars and zero the other ones
    optional_mask :=  (c, n, opts) -> When(IsBound(opts.trueSVSemantics) and opts.trueSVSemantics and not(n=2),
        vuperm_2x64f_spu(c, List([1..2],x->When(x<=n, x, 128))),
        c),

    # load contiguous with unaligned loads
    loadc := (self, sv, opts) >> Cond(sv=1,
        (y,x) -> Checked(ObjId(x)=nth, assign(y, vuperm_2x64f_spu(promote_spu2x64f(x, 0), [1, 128]))),
        sv=2,
        (y,x) -> Checked(ObjId(x)=nth, assign(y, vperm_2x64f_spu(promote_spu2x64f(x, 0), promote_spu2x64f(nth(x.loc, x.idx+1), 0), [1, 3]))),
        Error("bad arguments for loadc")),

    #load contiguous + known alignment -> using 2 aligned load to be smarter
    loadc_align := (self, sv, align, opts)>>
    ((y,x,addr) -> let(
        v1 := nth(nth(x,add(addr,-align)).toPtr(TVect(TReal, 2)),0),
        v2 := nth(nth(x,add(addr,2-align)).toPtr(TVect(TReal, 2)), 0),
        m := x-> self.optional_mask(x, sv, opts),
        Cond(align=0,
            assign(y, m(v1)),
            assign(y, vperm_2x64f_spu(v1, v2, List([1..2], x-> When(x<=sv, align+x, 128))))
            ))),

###########################################################################
###### </Vas>
###########################################################################

    svstore := [[ # store using subvectors of length 1
               (y,x,opts) ->        assign(y[1], extract_spu2x64f(x, 0)),

               (y,x,opts) -> chain( assign(y[1], extract_spu2x64f(x, 0)),
                               assign(y[2], extract_spu2x64f(x, 1)) ),
                ],
                [ # store using subvectors of length 2
                (y,x,opts) -> chain( assign(y[1],                      extract_spu2x64f(x, 0)),
                                assign(nth(y[1].loc, y[1].idx+1), extract_spu2x64f(x, 1)) ),
               ]],

    storec := [    # store contiguous unaligned
              (y,x) ->        assign(y, extract_spu2x64f(x, 0)),
              (y,x) -> chain( assign(y,                   extract_spu2x64f(x, 0)),
                              assign(nth(y.loc, y.idx+1), extract_spu2x64f(x, 1)) ),
              ],
    reverse := (y,x) -> assign(vref(y,0,2), vuperm_2x64f_spu(vref(x,0,2), [2, 1])),

    # Complex multiplication
    mul_cx := (self, opts) >> Cond(
               opts.vector.SIMD in ["SSE"],
               Error("opts.vector.SIMD is in SSE, while I am in SPU. You've got problems!"),
               (y, x, c) ->
                    let(u1 := var.fresh_t("U", TVectDouble(2)), u2 := var.fresh_t("U", TVectDouble(2)),
                        u3 := var.fresh_t("U", TVectDouble(2)), u4 := var.fresh_t("U", TVectDouble(2)),
                        decl([u1, u2, u3, u4], chain(
                            assign(u1, mul(x, vuperm_2x64f_spu(c, [1,1]))),
                            assign(u2, vuperm_2x64f_spu(x, [2,1])),
                            assign(u3, mul(u2, vuperm_2x64f_spu(c, [2,2]))),
                            assign(u4, vuperm_2x64f_spu(u3, [2,1])),
                            assign(y, add(u1, u4)))))
                ),

    #----------------------------------------
    # support for VS and VS.transpose()

    shl1 :=  (y,x,opts) -> assign(y, vperm_2x64f_spu(x[1], x[2], [2, 3])),
    shl2 :=  (y,x,opts) -> assign(y, vperm_2x64f_spu(x[1], x[2], [2, 128])),
    shr1 :=  (y,x,opts) -> assign(y, vperm_2x64f_spu(x[1], x[2], [4, 1])),
    shr2 :=  (y,x,opts) -> assign(y, vperm_2x64f_spu(x[1], x[2], [128, 1])),

    # support for VO1dsJ(n, v)
    shrev := (y,x,opts) -> assign(y, vperm_2x64f_spu(x[1], x[2], [1, 2])),

    swap_cx := (y, x, opts) -> assign(y, vuperm_2x64f_spu(x, [2,1])),
    RCVIxJ2 := (y, x, opts) -> assign(y, vuperm_2x64f_spu(x, [2,1]))

));

Class(ppu_4x32f, spu_4x32f, rec(
    includes := () -> ["<omega32.h>", "<vec_types.h>", "<altivec.h>", "<spu2vmx.h>"],
    backendConfig := rec(
                        profile := default_profiles.linux_cellPPU_gcc,
                        measureFunction := _StandardMeasureVerify
                     )
));

Class(spu_8x16i, spu_common, rec(
    active          := false,
    isFixedPoint    := true,
    saturatedArithmetic := true,
    info            := "SPU 8 x 16-bit int",
    v               := 8,
    ctype           := "short int",
    vtype           := "vector signed short",
    stype           := "__attribute__ ((aligned(16))) short int",
    # NOTE: What is svtype?
    svtype          := ["short int", "int", "signed long long", "signed long long"],
    instr           := [ vperm_8x16i_spu, vuperm_8x16i_spu ],
    bits            := 16,
    isFloat         := false,
    isFix           := true,
    includes        := () -> ["<include/omega16i.h>", "<include/mm_malloc.h>", "<spu_intrinsics.h>"],
    infix_op        := false,
    infix_assign    := true,
    vzero           := vzero_8x16i,
    vconstv         := "(vector signed short)",
    vconst1         := "(vector signed short)spu_splats",
    #NOTE: What is vconsthex? - it's to set 8 values. Value in unparser should take care of it?
    #vconsthex       := p->Chain(Print("_mm_set_epi16("), PrintCS(p), Print(")")),
    #NOTE: what are the following two?
    #vconst1init     := c -> PrintCS(List([1..8], i->c)),
    #vconst          := vconstprfp,
    splopts         := rec(customDataType := "signed short"),

    #NOTE: These REALLY need to be rewritten. They're correct, but end up doing
    #extra loads, shifts, rotates like nobody's business.
#     loadc := [      # load contiguous unaligned
#        (y,x) -> Checked(ObjId(x)=nth, assign(y, vuperm_8x16i_spu(promote_spu8x16i(x, 0), [1, 128, 128, 128, 128, 128, 128, 128]))),
#        (y,x) -> Checked(ObjId(x)=nth, assign(y, vperm_8x16i_spu(promote_spu8x16i(x, 0), promote_spu8x16i(nth(x.loc, x.idx+1), 0), [1,5,128,128,128,128,128,128]))),
#        (y,x) -> Checked(ObjId(x)=nth, assign(y, vuperm_8x16i_spu(vloadu8_spu8x16i(x), [1,2,3,128,128,128,128,128]))),
#        (y,x) -> Checked(ObjId(x)=nth, assign(y, vuperm_8x16i_spu(vloadu8_spu8x16i(x), [1,2,3,  4,128,128,128,128]))),
#        (y,x) -> Checked(ObjId(x)=nth, assign(y, vuperm_8x16i_spu(vloadu8_spu8x16i(x), [1,2,3,  4,  5,128,128,128]))),
#        (y,x) -> Checked(ObjId(x)=nth, assign(y, vuperm_8x16i_spu(vloadu8_spu8x16i(x), [1,2,3,  4,  5,  6,128,128]))),
#        (y,x) -> Checked(ObjId(x)=nth, assign(y, vuperm_8x16i_spu(vloadu8_spu8x16i(x), [1,2,3,  4,  5,  6,  7,128]))),
#        (y,x) -> Checked(ObjId(x)=nth, assign(y, vloadu8_spu8x16i(x)))
#     ],

###########################################################################
###### Vas, please check
###########################################################################

    # keep the n lower scalars and zero the other ones
    optional_mask :=  (c, n, opts) -> When(IsBound(opts.trueSVSemantics) and opts.trueSVSemantics and not(n=8),
        vuperm_8x16i_spu(c, List([1..8],x->When(x<=n, x, 128))),
        c),

    #NOTE: These REALLY need to be rewritten. They're correct, but end up doing
    #extra loads, shifts, rotates like nobody's business.
    # load contiguous with unaligned loads
    loadc := (self, sv, opts) >> ((y,x) -> assign(y, self.optional_mask(vloadu8_spu8x16i(x.toPtr(TVect(TReal, 8))), sv, opts))),

    #load contiguous + known alignment -> using 2 aligned load to be smarter
    loadc_align := (self, sv, align, opts)>>
    ((y,x,addr) -> let(
        v1 := nth(nth(x,add(addr,-align)).toPtr(TVect(TReal, 8)),0),
        v2 := nth(nth(x,add(addr,8-align)).toPtr(TVect(TReal, 8)), 0),
        m := x-> self.optional_mask(x, sv, opts),
        Cond(align=0,
            assign(y, m(v1)),
            assign(y, vperm_8x16i_spu(v1, v2, List([1..8], x-> When(x<=sv, align+x, 128))))
            ))),

###########################################################################
###### </Vas>
###########################################################################

    #NOTE: See if there's a better way to do this. Also look at assembly to see how good/bad this is.
    storec := [    # store contiguous unaligned
            (y,x) ->        assign(y, extract_spu8x16i(x, 0)),

            (y,x) -> chain( assign(y,                   extract_spu8x16i(x, 0)), assign(nth(y.loc, y.idx+1), extract_spu8x16i(x, 1))),

            (y,x) -> chain( assign(y,                   extract_spu8x16i(x, 0)), assign(nth(y.loc, y.idx+1), extract_spu8x16i(x, 1)),
                            assign(nth(y.loc, y.idx+2), extract_spu8x16i(x, 2))),

            (y,x) -> chain( assign(y,                   extract_spu8x16i(x, 0)), assign(nth(y.loc, y.idx+1), extract_spu8x16i(x, 1)),
                            assign(nth(y.loc, y.idx+2), extract_spu8x16i(x, 2)), assign(nth(y.loc, y.idx+3), extract_spu8x16i(x, 3))),

            (y,x) -> chain( assign(y,                   extract_spu8x16i(x, 0)), assign(nth(y.loc, y.idx+1), extract_spu8x16i(x, 1)),
                            assign(nth(y.loc, y.idx+2), extract_spu8x16i(x, 2)), assign(nth(y.loc, y.idx+3), extract_spu8x16i(x, 3)),
                            assign(nth(y.loc, y.idx+4), extract_spu8x16i(x, 2))),

            (y,x) -> chain( assign(y,                   extract_spu8x16i(x, 0)), assign(nth(y.loc, y.idx+1), extract_spu8x16i(x, 1)),
                            assign(nth(y.loc, y.idx+2), extract_spu8x16i(x, 2)), assign(nth(y.loc, y.idx+3), extract_spu8x16i(x, 3)),
                            assign(nth(y.loc, y.idx+4), extract_spu8x16i(x, 2)), assign(nth(y.loc, y.idx+5), extract_spu8x16i(x, 5))),

            (y,x) -> chain( assign(y,                   extract_spu8x16i(x, 0)), assign(nth(y.loc, y.idx+1), extract_spu8x16i(x, 1)),
                            assign(nth(y.loc, y.idx+2), extract_spu8x16i(x, 2)), assign(nth(y.loc, y.idx+3), extract_spu8x16i(x, 3)),
                            assign(nth(y.loc, y.idx+4), extract_spu8x16i(x, 2)), assign(nth(y.loc, y.idx+5), extract_spu8x16i(x, 5)),
                            assign(nth(y.loc, y.idx+6), extract_spu8x16i(x, 2))),

            (y,x) -> chain( assign(y,                   extract_spu8x16i(x, 0)), assign(nth(y.loc, y.idx+1), extract_spu8x16i(x, 1)),
                            assign(nth(y.loc, y.idx+2), extract_spu8x16i(x, 2)), assign(nth(y.loc, y.idx+3), extract_spu8x16i(x, 3)),
                            assign(nth(y.loc, y.idx+4), extract_spu8x16i(x, 2)), assign(nth(y.loc, y.idx+5), extract_spu8x16i(x, 5)),
                            assign(nth(y.loc, y.idx+6), extract_spu8x16i(x, 2)), assign(nth(y.loc, y.idx+7), extract_spu8x16i(x, 5))),
        ],
    reverse := (y,x) -> assign(vref(y,0,8), vuperm_8x16i_spu(vref(x,0,4), [8, 7, 6, 5, 4, 3, 2, 1]))
));

SIMD_ISA_DB.addISA(spu_4x32f);
SIMD_ISA_DB.addISA(spu_2x64f);
SIMD_ISA_DB.addISA(ppu_4x32f);


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

_VConstructR := [VBase, VRCDiag, BlockVPerm, VPerm, VTensor, VDiag_x_I, VDiag, VGath, VScat, VScatAcc, VPrm_x_I, VGath_sv, VBlk, VContainer, RCVDiag, VO1dsJ, VScale, VGath_zero, VReplicate, VIxJ2, VJxI];
_VConstructL := [VBase, VRCDiag, BlockVPerm, VPerm, VTensor, VDiag_x_I, VDiag, VGath, VScat, VScatAcc, VPrm_x_I, VScat_sv, VScat_svAcc, VScat_pc, VBlk, VContainer, RCVDiag, VReplicate, VO1dsJ, VScale, VGath_zero, VIxJ2, VJxI];

_VTags := [AVecReg, AVecRegCx];

_v_divides_cols := x -> _divides(getV(x), Cols(x));
_v_divides_rows := x -> _divides(getV(x), Rows(x));
_v_tagged_ntR := x->IsNonTerminal(x) and x.hasAnyTag(_VTags) and _divides(x.getAnyTag(_VTags).v, Rows(x));
_v_tagged_ntL := x->IsNonTerminal(x) and x.hasAnyTag(_VTags) and _divides(x.getAnyTag(_VTags).v, Cols(x));

Scat.svVariant := VScat_sv;
ScatAcc.svVariant := VScat_svAcc;
Gath.svVariant := VGath_sv;

Class(RulesVec, RuleSet);
RewriteRules(RulesVec, rec(
    ############################################################
    # Below 2 rules are required by SAR and autolib
    # NOTE: They will cause problems when mixing scalar and vector code,
    # because the entire code (including the scalar portion) will reside inside
    # VContainer, and gath/scat in scalar code will be converted to vgath/vscat
    #
    VContainer__Gath_Scat := Rule(@@(1, [Gath, Scat, ScatAcc], (e,cx)->
        cx.isInside(VContainer) and IsSIMD_ISA(Last(cx.VContainer).isa)
        and let(sg := e, v := Last(cx.VContainer).isa.getV(), vsg := sg.svVariant(sg.func, v, 1), sg.dims() = vsg.dims())
        and (not IsBound(cx.VTensor) or cx.VTensor = [])
        and (not IsBound(cx.VJamL) or cx.VJamL = [])
        and (not IsBound(cx.VJamR) or cx.VJamR = [])
        and (not IsBound(cx.VJam1) or cx.VJam1 = [])
        and (not IsBound(cx.VJam ) or cx.VJam  = [])
        # FF, NOTE: this is a terrible hack to solve a contradiction: How to prevent the rader perm
	#            from turning into VGath. Need to find a better solution.
        and (not ObjId(e.func) in [RM, RR])),
     (e,cx) -> let(
	 #inside_vrc := (IsBound(cx.vRC) and cx.vRC <> []) or
	 #              (IsBound(cx.opts.assumeOuter_vRC) and cx.opts.assumeOuter_vRC),
	 v := Last(cx.VContainer).isa.getV(), # / Cond(inside_vrc, 2, 1),
	 # YSV: isa.getV() now takes care of 'inside_vrc', which is thus commented out,
	 #      and can be removed.
         e.svVariant(e.func, v, 1))),

    VContainer__Diag := Rule(@@(1, [Diag], (e,cx) ->
        cx.isInside(VContainer) and IsSIMD_ISA(Last(cx.VContainer).isa)
        and not cx.isInside(VTensor)
        and not cx.isInside("VJamL")
        and not cx.isInside("VJamR")
        and not cx.isInside("VJam1")
        and not cx.isInside("VJam")
	and let(
	    #inside_vrc := cx.isInside(vRC) or (IsBound(cx.opts.assumeOuter_vRC) and cx.opts.assumeOuter_vRC),
	    v := Last(cx.VContainer).isa.getV(), # / Cond(inside_vrc, 2, 1),
	    Rows(e) mod v = 0)),
     (e,cx) -> let(
	 #inside_vrc := (IsBound(cx.vRC) and cx.vRC <> []) or
	 #              (IsBound(cx.opts.assumeOuter_vRC) and cx.opts.assumeOuter_vRC),
	 v := Last(cx.VContainer).isa.getV(), # / Cond(inside_vrc, 2, 1),
	 # YSV: isa.getV() now takes care of 'inside_vrc', which is thus commented out,
	 #      and can be removed.
         VDiag(e.element, v))),

# YSV: if VContainer is used, then I need the following (in particular the VDiag)
    VContainer_XXX := ARule(Compose, [@(1, VContainer), @(2, [VDiag, Diag, Gath, Scat, ScatAcc])],
        e -> [ CopyFields(@(1).val, rec(_children := [@(1).val._children[1] * @(2).val])) ]),

    XXX_VContainer := ARule(Compose, [@(1, [VDiag, Diag, Gath, Scat, ScatAcc]), @(2, VContainer)],
        e -> [ CopyFields(@(2).val, rec(_children := [@(1).val * @(2).val._children[1]])) ]),

    VContainer_VContainer := Rule( [@(1, VContainer), @(2, VContainer, x -> x.isa = @(1).val.isa and x.isa.isCplx()=@(1).val.isa.isCplx())],
        e -> e.child(1)),

#    Gath_Scat_ScatGath_XXX := ARule(Compose, [ @(1, Scat), @(2, ScatGath) ],
##    Gath_Scat_ScatGath_XXX := ARule(Compose, [@@(1, [Scat, Gath], (e,cx)->IsBound(cx.VContainer) and cx.VContainer <> []), @(2, ScatGath)],
#        (e, cx)->let(sg := @@(1).val, v := Last(cx.VContainer).v,
#            Concat(Error(Caught), [ sg.svVariant(sg.func, vrc.v, 1)]) )),

#############################
    # Composition of subvector and full vector accesses
    # The result is VScat_sv that has rem=0 guaranteed, i.e., does not use partially filled vectors
    VScat_sv__VScat := ARule(Compose, [@(1,[VScat_sv, VScat_svAcc]), @(2, [VScat, VScatAcc])],
    e -> let(func := @(1).val.func,  vfunc := @(2).val.func,
             v    := getV(@(1).val), sv    := @(1).val.sv,
	     oid  := Cond(ObjId(@(1).val)=VScat_svAcc or ObjId(@(2).val)=VScatAcc, VScat_svAcc, VScat_sv),
        [ oid(fCompose(func, fTensor(vfunc, fId(v/sv))), v, sv, 0) ])),

    # The result is VGath_sv that has rem=0 guaranteed, i.e., does not use partially filled vectors
    VGath__VGath_sv := ARule(Compose, [@(1,VGath), @(2,VGath_sv)],
    e -> let(func := @(2).val.func,  vfunc := @(1).val.func,
             v    := getV(@(2).val), sv    := @(2).val.sv,
         [ VGath_sv(fCompose(func, fTensor(vfunc, fId(v/sv))), v, sv, 0) ])),

    # Composition of subvector and scalar accesses
    Scat__VScat_sv := ARule(Compose, [@(1, Scat), @(2, VScat_sv)], e -> let(s := @(2).val,
        [ VScat_sv(fCompose(@(1).val.func, fTensor(s.func, fId(s.sv))), s.v, 1, s.rem) ])),

    ScatAcc__VScat_sv := ARule(Compose, [@(1, ScatAcc), @(2, VScat_sv)], e -> let(s := @(2).val,
        [ VScat_svAcc(fCompose(@(1).val.func, fTensor(s.func, fId(s.sv))), s.v, 1, s.rem) ])),

    VGath_sv__Gath := ARule(Compose, [@(1, VGath_sv), @(2, Gath)], e -> let(s := @(1).val,
        [ VGath_sv(fCompose(@(2).val.func, fTensor(s.func, fId(s.sv))), s.v, 1, s.rem) ])),

    # Subvector access --> full vector access
    VScat_sv_to_VScat := Rule(@(1,VScat_sv,e->getV(e)=e.sv), e->VScat(@(1).val.func, @(1).val.v)),
    VGath_sv_to_VGath := Rule(@(1,VGath_sv,e->getV(e)=e.sv), e->VGath(@(1).val.func, @(1).val.v)),
    VScat_svAcc_to_VScatAcc := Rule(@(1,VScat_svAcc,e->getV(e)=e.sv), e->VScatAcc(@(1).val.func, @(1).val.v)),

    # Increase granularity if possible (when function is fId(n) or fTensor(X, fId(n)))
    GathScat_sv_fTensor := Rule([@(1, [VGath_sv, VScat_sv, VScat_svAcc]),
         [@(2,fTensor), ..., [fId, @(3).cond(e -> Gcd(@(1).val.v/@(1).val.sv, EvalScalar(e)) > 1 )]]],
     e -> let(v := @(1).val.v, sv := @(1).val.sv,
          n := EvalScalar(@(3).val),   gcd := Gcd(v / sv, n),
	  rem := Cond(@(1).val.rem _is Unk, Unk(TInt), @(1).val.rem/gcd),
          ObjId(@(1).val)(fTensor(DropLast(@(2).val.children(), 1), fId(n/gcd)), v, sv*gcd, rem))),

    GathScat_sv_HofTensor := Rule([@(1, [VGath_sv, VScat_sv, VScat_svAcc]),
                        [fCompose, @(9, H, e->_divides(@(1).val.sv, e.params[3]) and e.params[4]=1),
			    [@(2,fTensor), ..., [fId, @(3).cond(e -> Gcd(@(1).val.v/@(1).val.sv, EvalScalar(e)) > 1)]]]],
     e -> let(v := @(1).val.v, sv := @(1).val.sv,
          n := EvalScalar(@(3).val),   gcd := Gcd(v / sv, n), hp := @(9).val.params,
	  rem := Cond(@(1).val.rem _is Unk, Unk(TInt), @(1).val.rem/gcd),
          ObjId(@(1).val)(fCompose(H(hp[1]/gcd, hp[2]/gcd, hp[3]/gcd, 1),  fTensor(DropLast(@(2).val.children(), 1), fId(n/gcd))), v, sv*gcd, rem))),

    GathScat_sv_fId := Rule([@(1, [VGath_sv, VScat_sv, VScat_svAcc]),
	                         [fId, @(3).cond(e -> Gcd(@(1).val.v/@(1).val.sv, EvalScalar(e)) = @(1).val.v)]],
     e -> let(v := @(1).val.v, sv := @(1).val.sv,
          n := EvalScalar(@(3).val),   gcd := Gcd(v / sv, n),
	  rem := Cond(@(1).val.rem _is Unk, Unk(TInt), @(1).val.rem/gcd),
          ObjId(@(1).val)(fId(n/gcd), v, sv*gcd))),

    # Convert objects that are adjacent to vector constructs to vector objects
    #

    VConstruct_Gath := ARule(Compose, [ @(1, _VConstructL), @(2, Gath, e -> _v_divides_cols(@(1).val)) ],
        e -> [ @(1).val, VGath_sv(@(2).val.func, getV(@(1).val), 1) ]),
    VConstruct_Scat := ARule(Compose, [ @(1, _VConstructL), @(2, Scat, e -> _v_divides_cols(@(1).val)) ],
        e -> [ @(1).val, VScat_sv(@(2).val.func, getV(@(1).val), 1) ]),
    VConstruct_ScatAcc := ARule(Compose, [ @(1, _VConstructL), @(2, ScatAcc, e -> _v_divides_cols(@(1).val)) ],
        e -> [ @(1).val, VScat_svAcc(@(2).val.func, getV(@(1).val), 1) ]),
    VConstruct_Diag := ARule(Compose, [ @(1, _VConstructL), [Diag, @(2).cond(e -> _v_divides_cols(@(1).val))] ],
        e -> [ @(1).val, VDiag(@(2).val, getV(@(1).val)) ]),

    Scat_VConstruct := ARule(Compose, [ @(1, Scat), @(2, _VConstructR, _v_divides_rows) ],
        e -> [ VScat_sv(@(1).val.func, getV(@(2).val), 1), @(2).val ]),
    ScatAcc_VConstruct := ARule(Compose, [ @(1, ScatAcc), @(2, _VConstructR, _v_divides_rows) ],
        e -> [ VScat_svAcc(@(1).val.func, getV(@(2).val), 1), @(2).val ]),
    Gath_VConstruct := ARule(Compose, [ @(1, Gath), @(2, _VConstructR, _v_divides_rows) ],
        e -> [ VGath_sv(@(1).val.func, getV(@(2).val), 1), @(2).val ]),
    Diag_VConstruct := ARule(Compose, [ [Diag, @(1)], @(2, _VConstructR, _v_divides_rows) ],
        e -> [ VDiag(@(1).val, getV(@(2).val)), @(2).val ]),

    VGath_Prm := ARule(Compose, [ @(1, VGath), @(2, Prm) ],
        e -> [ @(1).val, VGath_sv(@(2).val.func, getV(@(1).val), 1) ]),
    Prm_VScat := ARule(Compose, [ @(1, Prm), @(2, VScat) ],
        e -> [ VScat_sv(@(1).val.func.transpose(), getV(@(2).val), 1), @(2).val ]),
    Prm_VScatAcc := ARule(Compose, [ @(1, Prm), @(2, VScatAcc) ],
        e -> [ VScat_svAcc(@(1).val.func.transpose(), getV(@(2).val), 1), @(2).val ]),

    # NOTE: RowVec must be of constant size, fix this limitation
    VRowVec := Rule([@(1, VTensor), @(2, RowVec, e -> not IsSymbolic(e.element.domain()))],
        e -> VTensor(@(2).val.toDiagBlk(), @(1).val.vlen)),

    VColVec := Rule([@(1, VTensor), @(2, ColVec, e -> not IsSymbolic(e.element.domain()))],
        e -> VTensor(@(2).val.toDiagBlk(), @(1).val.vlen)),

    VPrm_x_I := Rule([@(1, VTensor), @(2, Prm)], e->VPrm_x_I(@(2).val.func, @(1).val.vlen)),
    VDiag_x_I := Rule([@(1, VTensor), @(2, Diag)], e->VDiag_x_I(@(2).val.element, @(1).val.vlen)),
    RCDiag_x_I := Rule([@(1, VTensor), @(2, RCDiag)], e->VRCDiag(VDup(@(2).val.element, @(1).val.vlen), @(1).val.vlen)),

    VPrm_x_I_Id := ARule(Compose,[@(1,VPrm_x_I), [Prm, fId]],e->[@(1).val]),
    Id_VPrm_x_I := ARule(Compose, [[Prm, fId], @(1,VPrm_x_I)], e->[@(1).val]),

    # Remove identity gathers / scatters
    RemGSP := Rule([@(1, [Gath, Scat, Prm, VGath, VScat]), @(2, fId)], e->I(@(1).val.dims()[1])),
    RemIL  := ARule(Compose, [ @(1), @(2, I) ], e -> [ @(1).val ]),
    RemIR  := ARule(Compose, [ @(1, I), @(2) ], e -> [ @(2).val ]),

    #   H rules
    VGath_sv_H := Rule( [@(1,VGath_sv), [H,
            @(2).cond(e->not IsSymbolic(e)),
            @(3).cond(e->_divides(getV(@(1).val), e)),
            @(4).cond(e->_divides(getV(@(1).val), e)), _1]],
	e -> let(
            d := e.v / e.sv,
	    rmod := EvalScalar(@(2).val) mod getV(@(1).val),
	    Cond(rmod=0,
		 VGath(H(@(2).val/d, @(3).val/d, @(4).val/d, 1), e.v),
		 VGath_pc(@(2).val*e.sv, @(3).val*e.sv, @(4).val*e.sv, e.v))
	)),

    VScat_sv_H := Rule( [@(1, [VScat_sv, VScat_svAcc]), [H,
            @(2).cond(e->not IsSymbolic(e)),
            @(3).cond(e->_divides(getV(@(1).val), e)),
            @(4).cond(e->_divides(getV(@(1).val), e)), _1]],
	e -> let(
            d := e.v / e.sv,
	    oid   := Cond(ObjId(@(1).val)=VScat_sv, VScat, VScatAcc),
	    oidpc := Cond(ObjId(@(1).val)=VScat_sv, VScat_pc, VScat_pcAcc),
	    rmod := EvalScalar(@(2).val) mod getV(@(1).val),
	    Cond(rmod=0,
		 oid(H(@(2).val/d, @(3).val/d, @(4).val/d, 1), e.v),
		 oidpc(@(2).val*e.sv, @(3).val*e.sv, @(4).val*e.sv, e.v))
	)),

   Scat_VScat_pc_to_VScat_sv := ARule(Compose, [@(1, Scat), @(2, VScat_pc, x -> x.N=x.n and x.ofs = 0)],
       e -> [VScat_sv(@(1).val.func, @(2).val.v, 1)]),

   VGath_pc_Gath_H := ARule(Compose, [@(1, VGath_pc), [Gath, [@(2, H), @, @, @, _1]]],
       e -> [VGath_pc(@(2).val.params[1], @(1).val.n, @(2).val.params[3], @(1).val.v)] ),

   Scat_H_VScat_pc := ARule(Compose, [[Scat, [@(1, H), @, @, @, _1]], @(2, [VScat_pc, VScat_pcAcc])],
       e -> [ObjId(@(2).val)(@(1).val.params[1], @(2).val.n, @(1).val.params[3], @(2).val.v)] ),

   ScatGath := Rule([@(1, VTensor), @(2, ScatGath)],
       e -> ScatGath(fTensor(@(2).val.sfunc, fId(@(1).val.vlen)), fTensor(@(2).val.gfunc, fId(@(1).val.vlen)))),

   Scat_Cvt := ARule(Compose, [@(1, Scat), @(2, Cvt, x -> IsSIMD_ISA(x.params[1].isa_to) and _divides(x.params[1].isa_to.v, Cols(@(1).val)) )],
      e -> let(isa := @(2).val.params[1].isa_to, [ VContainer(VScat_sv(@(1).val.func, isa.v, 1), isa), @(2).val ])),
   Cvt_Gath := ARule(Compose, [@(1, Cvt, x -> IsSIMD_ISA(x.params[1].isa_from)), @(2, Gath, x -> _divides(@(1).val.params[1].isa_from.v, Rows(x)))],
      e -> let(isa := @(1).val.params[1].isa_from, [ @(1).val, VContainer(VGath_sv(@(2).val.func, isa.v, 1), isa) ])),
));


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

Class(assign_cmd, ExpCommand, rec(
    op_in := self >> Set( List(Drop(self.args,1), ArgsExp)),
    op_out := self >> Set([self.args[1]]),
    op_inout := self >> Set([]),
    unparse := "assign_cmd",
));

Class(neg_cmd, assign_cmd);

Class(fma_cmd,  assign_cmd, rec(exp_op := fma));
Class(fms_cmd,  assign_cmd, rec(exp_op := fms));
Class(nfma_cmd, assign_cmd, rec(exp_op := nfma));

# Examples:
#   assign(t1, add(t1, t2)) == assign_add(t1, t2),
#   assign(t1, add(t2, t3)) == chain(assign(t1, t2), assign_add(t2, t3)).
#

_threeop_bintab := WeakRef(tab());
_threeop_binlist := [];

ThreeOpFromBinOp := function(arg)
    local op;
    for op in arg do
        _threeop_bintab.(op.exp_op.__name__) := op;
        Add(_threeop_binlist, op.exp_op);
    od;
end;

ThreeOpFromBinOp(
   Class(add_cmd, assign_cmd, rec( exp_op := add )),
   Class(mul_cmd, assign_cmd, rec( exp_op := mul )),
   Class(sub_cmd, assign_cmd, rec( exp_op := sub ))
);

Class(ThreeOpRuleSet, RuleSet);
RewriteRules(ThreeOpRuleSet, rec(
    neg := Rule([assign, @(1), [neg, @(2)]], e -> neg_cmd(@(1).val, @(2).val)),

    fma  := Rule([assign, @(1), [ fma, @(2), @(3), @(4)]], e ->  fma_cmd(@(1).val, @(2).val, @(3).val, @(4).val)),
    fms  := Rule([assign, @(1), [ fms, @(2), @(3), @(4)]], e ->  fms_cmd(@(1).val, @(2).val, @(3).val, @(4).val)),
    nfma := Rule([assign, @(1), [nfma, @(2), @(3), @(4)]], e -> nfma_cmd(@(1).val, @(2).val, @(3).val, @(4).val)),

    binop := Rule([assign, @(1), [@(0,_threeop_binlist) , @(2), @(3)]], e -> _threeop_bintab.(@(0).val.__name__)(@(1).val, @(2).val, @(3).val)),
));

ThreeOpMacroUnparser_Mixin := rec(

    add_cmd := (self,o,i,is) >> Print(Blanks(i), self.prefixTTT("ADD", o.args[1].t, o.args[2].t, o.args[3].t, o.args), ";\n"),
    sub_cmd := (self,o,i,is) >> Print(Blanks(i), self.prefixTTT("SUB", o.args[1].t, o.args[2].t, o.args[3].t, o.args), ";\n"),
    mul_cmd := (self,o,i,is) >> Print(Blanks(i), self.prefixTTT("MUL", o.args[1].t, o.args[2].t, o.args[3].t, o.args), ";\n"),
    neg_cmd := (self,o,i,is) >> Print(Blanks(i), self.prefixTT( "NEG", o.args[1].t, o.args[2].t, o.args), ";\n"),

    fma_cmd  := (self,o,i,is) >> Print(Blanks(i), self.prefixTTTT("FMA", o.args[1].t, o.args[2].t, o.args[3].t, o.args[4].t, o.args), ";\n"),
    fms_cmd  := (self,o,i,is) >> Print(Blanks(i), self.prefixTTTT("FMS", o.args[1].t, o.args[2].t, o.args[3].t, o.args[4].t, o.args), ";\n"),
    nfma_cmd := (self,o,i,is) >> Print(Blanks(i), self.prefixTTTT("NFMA", o.args[1].t, o.args[2].t, o.args[3].t, o.args[4].t, o.args), ";\n"),

);

DoThreeOp := function(c, opts)
    c := MarkDefUse(c);
    c := BinSplit(c);
    c := ThreeOpRuleSet(c);
    return c;
end;

ThreeOp_CS := Concatenation(BaseIndicesCS, [
    MarkDefUse, #
    (c, opts) -> CopyPropagate(c, CopyFields(opts, rec(autoinline := true))),
    MarkDefUse, #
    (c, opts) -> CopyPropagate(c, CopyFields(opts, rec(autoinline := true))),
    DoThreeOp,
    Compile.declareVars
]);

Class(ThreeOpUnparser, CMacroUnparserProg, rec(
    add_cmd := (self,o,i,is) >> self.prefixTTT("ADD", o.args[1].t, o.args[2].t, o.args[3], o.args),

    sub_cmd := (self,o,i,is) >> self.prefixTTT("SUB", o.args[1].t, o.args[2].t, o.args[3], o.args),

    neg_cmd := (self,o,i,is) >> self.prefixT("NEG", o.t, o.args),

    # first argument of mul_cmd is the result, 2nd and 3rd are operands
    mul_cmd := (self,o,i,is) >> When(Length(o.args)=3,
            let(
            # check if constant ended up in slot #3
            a := When(IsValue(o.args[3]), o.args[3], o.args[2]),
            b := When(IsValue(o.args[3]), o.args[2], o.args[3]),
            When(not (IsValue(a) or (IsVar(a) and IsBound(a.value))),
                # <a> is not a constant
                self.prefixTTT("MUL", o.args[1], a, b, [o.args[1],a,b]),
                #self.prefix(Concat("MUL_", self._pfx(a.t), "_", self._pfx(b.t)), [o.args[1],a,b]),
                # <a> is a constant
                let(fmt := self._const( When(IsValue(a), a, a.value) ),
                    self.prefix(Concat("MUL_", o.args[1]._pfx(o.args[1].t), "CNST", "_", self._pfx(b.t)), [o.args[1], GetExponent(a.value), GetMantissa(a.value), b]))))),
#                   self.prefix(Concat("MUL_", o.args[1]._pfx(o.args[1].t), fmt[1], "_", self._pfx(b.t)), [o.args[1], a, b]))))),
#                        # check if constant is special, and does not go into MUL args
#                        # for example fmt[1]="I" denotes sqrt(-1), one such constant
#                        #When(fmt[2]=[], [b], [a, b])))))),

#    add_cmd := (self, o, i, is) >> Print(Blanks(i), self.prefix(self, i, is), ";\n"),
#    sub_cmd := ~.add_cmd

    prefixT := (self, funcname, t, args) >>
        self.prefix(Concat(funcname, "_", self._pfx(t)), args),

    prefixTT := (self, funcname, t1, t2, args) >>
        self.prefix(Concat(funcname, "_", self._pfx(t1), "_", self._pfx(t2)), args),

    prefixTTT := (self, funcname, t1, t2, t3, args) >>
        self.prefix(Concat(funcname, "_", self._pfx(t1), "_", self._pfx(t2), "_", self._pfx(t3)), args),

    _pfx := (self, t) >> Cond(
        ObjId(t) = T_Complex, "CPX",
        ObjId(t) = T_Real,    "FLT",
        t = TComplex, "CPX",
        t = TReal,    "FLT",
        t = TInt or ObjId(t) in [TArray, TPtr], "INT",
        t = TUnknown, "UNK",
        ObjId(t) = TSym, "SYM",
        IsVecT(t), Cond(
            t.t = TReal, Concat("FV",StringInt(t.size)),
            t.t = TComplex, Concat("FC",StringInt(t.size)),
            t.t = TInt, Concat("IV",StringInt(t.size)),
            Error("Can't handle type ", t)
        ),
        Error("Can't handle type ", t)),

    # returns a tuple [suffix, args], where suffix is used for MUL_XXX or C_XXX,
    # and args are additional parameters into C_XXX
    _const := (self,o) >> Cond(
        o.t = TReal and IsCyc(o.v),        ["FLT", [ReComplex(Complex(o.v))]],
        o.t = TReal,                       ["FLT", [o.v]],
        o.t = TInt,                        ["INT", [o.v]],
        o.t = TUnknown,                    ["INT", [o.v]], # NOTE: there is a bug that creates V(0) with TUnknown
        o.t = TString,                     ["STR", [o.v]],
        o.t = TBool,                       ["INT", [When(o.v in [true, 1], 1, 0)]],
        Error("Don't know how to handle constant of type ", o.t)
    ),
));

EncodeFloat := function(f)
    local a, exp, mant, tmp, flag;

    #shift by max exponent, which is +/- 127 in IEEE floating point
    if (IntDouble(f)=0) then
        exp := Log2Int(IntDouble(f*2^127))-127;
    else
        exp := Log2Int(IntDouble(f));
    fi;

    tmp := f/(2^exp);
    mant := tmp - IntDouble(tmp);
    mant := IntDouble(mant * 2^23);

        a:= rec(exp:= exp+127, mant:=mant);

    return a;
end;

DecodeToFloat := function(exp, mant)
    local f;

    if (mant<0) then
        mant := -mant;
        f := mant / (2^23);
        f := f + 1;
        f:= f * (2^(exp-127));
        f:= -f;
    else
        f := mant / (2^23);
        f := f + 1;
        f:= f * (2^(exp-127));
    fi;

    return f;
end;

GetExponent := function(f)
    local exp;

    #shift by max exponent, which is +/- 127 in IEEE floating point
    if (IntDouble(f)=0) then
        exp := Log2Int(IntDouble(f*2^127))-127;
    else
        exp := Log2Int(IntDouble(f));
    fi;

    return exp + 127;
end;

GetMantissa := function(f)
    local a, exp, mant, tmp, flag;

    #shift by max exponent, which is +/- 127 in IEEE floating point
    if (IntDouble(f)=0) then
        exp := Log2Int(IntDouble(f*2^127))-127;
    else
        exp := Log2Int(IntDouble(f));
    fi;

    tmp := f/(2^exp);
    mant := tmp - IntDouble(tmp);
    mant := IntDouble(mant * 2^23);

    return mant;
end;


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

Declare(OLCompose);

Class(OLCompose, Compose, rec(
    printSeparationChar := " o ",
    toOperator := self >> (vec -> Checked(Length(vec)=self.dims()[2], OLCompose(DropLast(self._children, 1)).toOperator()(Last(self._children).toOperator()(vec))))
));

Class(LinearSystem,  BaseMat, rec(
    abbrevs := [ (A, b) -> [A, b] ],
    new := (self, A, b) >> SPL(WithBases(self, rec(
        mat := A, rhs := b, dimensions := [Rows(A), Cols(A)]))),
    print := (self, i, is) >> Print("LinearSystem(", self.mat, ", ", self.rhs, ")"),
    dims  := self >> self.dimensions,
    rChildren := self >> [self.mat, self.rhs],
    rSetChild := rSetChildFields("mat", "rhs"),
    isReal := True
));

Class(Reduction,  BaseMat, rec(
    abbrevs := [ (N, op, idval, isSaturated) -> [N, op, idval, isSaturated] ],
    new := (self, N, op, idval, isSaturated) >> SPL(WithBases(self, rec(
        N := N, op := op, idval := idval, isSaturated := isSaturated, dimensions := [1, N]))),
    print := (self, i, is) >> Print("Reduction(", self.N, ", ", self.op, ", ", self.idval, ", ", self.isSaturated, ")"),
    dims  := self >> self.dimensions,
    rChildren := self >> [self.N, self.op, self.idval, self.isSaturated],
    rSetChild := rSetChildFields("N", "op", "idval", "isSaturated"),
    transpose := self >> self,
    isReal := True,
    toOperator := self >> (vec -> Checked(Length(vec)=self.dims()[2], [RulesStrengthReduce(FoldL(vec, self.op, self.idval)).v]))
));

Class(Induction,  BaseMat, rec(
    abbrevs := [ (N, op, initval) -> [N, op, initval] ],
    new := (self, N, op, initval) >> SPL(WithBases(self, rec(
        N := N, op := op, initval := initval, dimensions := [N, 1]))),
    print := (self, i, is) >> Print("Induction(", self.N, ", ", self.op, ", ", self.initval, ")"),
    dims  := self >> self.dimensions,
    rChildren := self >> [self.N, self.op, self.initval],
    rSetChild := rSetChildFields("N", "op", "initval"),
    isReal := True,
    recurrence := self >> ((j,v)->When(j=0, self.initval, self.op.at(self.recurrence()(j-1,v),v))),
    toOperator := self >> (vec -> Checked(Length(vec)=self.dims()[2], List([0..self.N-1], i->_unwrap(RulesStrengthReduce(self.recurrence()(i,vec[1]))))))
));

Class(PointWise, BaseMat, rec(
    abbrevs := [ (N, op) -> [N, op] ],
    new := (self, N, op) >> SPL(WithBases(self, rec(
        N := N, op := op, dimensions := [N, N]))),
    print := (self, i, is) >> Print("PointWise(", self.N, ", ", self.op, ")"),
    dims  := self >> self.dimensions,
    rChildren := self >> [self.N, self.op],
    rSetChild := rSetChildFields("N", "op"),
    transpose := self >> self,
    isReal := True,
    toOperator := self >> (vec -> Checked(Length(vec)=self.dims()[2], List([0..Length(vec)-1], j->RulesStrengthReduce(self.op.at(vec[j+1], j)).v)))
));

Class(BinOp, BaseMat, rec(
    abbrevs := [ (N, op) -> [N, op] ],
    new := (self, N, op) >> SPL(WithBases(self, rec(
        N := N, op := op, dimensions := [N, 2*N]))),
    print := (self, i, is) >> Print("BinOp(", self.N, ", ", self.op, ")"),
    dims  := self >> self.dimensions,
    rChildren := self >> [self.N, self.op],
    rSetChild := rSetChildFields("N", "op"),
    isReal := True,
    toOperator := self >> (vec -> Checked(Length(vec)=self.dims()[2], List([1..Length(vec)/2], j->RulesStrengthReduce(self.op.at(vec[j], vec[j+self.N])).v)))
));

Class(ScalarProd,  BaseMat, rec(
    abbrevs := [ (N, scalars) -> [N, scalars] ],
    new := (self, N, scalars) >> SPL(WithBases(self, rec(
        N := N, scalars := scalars, dimensions := [1, N]))),
    print := (self, i, is) >> Print("ScalarProd(", self.N, ", ", self.scalars, ")"),
    dims  := self >> self.dimensions,
    rChildren := self >> [self.N, self.scalars],
    rSetChild := rSetChildFields("N", "scalars"),
    isReal := True,
    toOperator := self >>(vec -> Checked(Length(vec)=self.dims()[2], [_unwrapVec(self.scalars) * vec]))
));

Class(Constant, BaseMat, rec(
    abbrevs := [ x -> [x]],
    new := (self, x) >> SPL(WithBases(self, rec(
    value := x, dimensions := [1,0]))),
    print := (self, i ,is) >> Print("Constant(",self.value,")"),
    dims := self >> self.dimensions,
    rChildren := self >> [self.value],
    rSetChild := rSetChildFields("value"),
    isReal := True
));

spiral.spl.IterVStack.toOperator :=
  self >> (vec -> Checked(Length(vec)=self.dims()[2],
                     Flat(List([0..self.var.range-1],
                        j->RulesStrengthReduce(
                               SubstVars(Copy(self._children[1]),
                                         rec((self.var.id):=V(j))
                               )
                           ).toOperator()(vec) )
                     )
                  )
          );

spiral.spl.HStack.toOperator :=
   self >> (vec -> Checked(Length(vec)=self.dims()[1]*Length(self.children()),
      	   	      Sum(List([1..Length(self.children())],
	                 i->self.children()[i].toOperator()(vec{[(i-1)*self.dims()[1] + 1..i*self.dims()[1]]})
		      ))
	 	   )
           );

spiral.spl.DirectSum.toOperator := self >> (vec -> Checked(Length(vec)=Sum(List(self.children(), i->i.dims()[2])),
    ApplyFunc(Concat, List([1..Length(self.children())],
        i->self.children()[i].toOperator()(vec{
            [Sum(List(self.children(){[1..i-1]}, j->j.dims()[2]))+1..Sum(List(self.children(){[1..i]}, j->j.dims()[2]))]
            })))));

spiral.spl.I.N := self >> self.params[1];

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: Read("file");

ext: gap

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


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

Class(loopc, loop_base, rec(
   __call__ := (self, loopvar, range, cnd, cmd) >> WithBases(self,
               rec(operations := CmdOps, cmd := cmd, var := loopvar, range := listRange(range), cnd := cnd)),
   rChildren := self >> [self.var, self.range, self.cnd, self.cmd],
   rSetChild := rSetChildFields("var", "range", "cnd", "cmd"),

   print := (self, i, si) >> Print(self.__name__, "(", self.var, ", ",
       self.range, ", ", self.cnd, ",\n", Blanks(i+si),
       self.cmd.print(i+si, si),
       Print("\n", Blanks(i), ")")),
   free := self >> self.cmd.free()
));

Class(creturn, throw);
Class(creturnCond, throw);

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

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


text:
signature Nuss(alphabet, answer) {

  answer nil(void);
  answer right(answer, alphabet);
  answer pair(alphabet, answer, alphabet);
  answer split(answer, answer);
  choice [answer] h([answer]);

}

algebra bpmax implements Nuss(alphabet = char, answer = int)
{
  int nil(void)
  {
    return 0;
  }

  int left(char c, int a)
  {
    return a;
  }

  int pair(char c, int m, char d)
  {
    return m + 1;
  }

  int split(int l, int r)
  {
    return l + r;
  }

  choice [int] h([int] l)
  {
    return list(maximum(l));
  }

}

grammar nussinov uses Nuss (axiom=start) {

  tabulated { start, bp }

  start = nil(EMPTY) |
          left(CHAR, start) |
          split(bp, start) # h ;

  bp = pair(CHAR, start, CHAR) with char_basepairing ;

}

instance bpmax = nussinov ( bpmax ) ;

max_stars_repo_name: adp-multi/adp-multi

max_issues_repo_name: adp-multi/adp-multi

max_forks_repo_name: adp-multi/adp-multi


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

_SKLR_ISA_TYPES := [
    [ T_Int(8),    "char"      ],
    [ T_Int(16),   "short"     ],
    [ T_Int(32),   "int"       ],
    [ T_Int(64),   "long long" ],
    [ T_UInt(8),   "unsigned char"      ],
    [ T_UInt(16),  "unsigned short"     ],
    [ T_UInt(32),  "unsigned int"       ],
    [ T_UInt(64),  "unsigned long long" ],
    [ T_Real(32),  "float"  ],
    [ T_Real(64),  "double" ],
    [ TInt,        "long"   ],
    [ TUInt,       "unsigned long"],
    [ TReal,       "double" ]
];

Class(SKLR, ISA, rec(
    __call__ := (self, el_t) >> WithBases(self, rec(
        t            := el_t,
	isSigned     := el_t.isSigned(),
	isFloat      := IsRealT(el_t),
	isFixedPoint := IsFixedPtT(el_t),
	operations   := ISAOps,
	print        := self >> Cond( self.t in List(_SKLR_ISA_TYPES, e -> e[1]),
	                                Print(self.id()), Print(self.__name__, "(", self.t, ")")),
	id           := self >> self.__name__ :: "_" :: el_t.strId(),
        ctype        := let( t := PositionProperty(_SKLR_ISA_TYPES, e -> e[1]=el_t),
                             Cond(t=false, "", _SKLR_ISA_TYPES[t][2])),
    )),

    isScalarISA  := true,

    v         := 1,
    getTags   := self >> [],
    getTagsCx := self >> [],
    toScalar  := self >> self,

    # this must be depricated:
    saturatedArithmetic := false,

));

IsSKLR_ISA := (x) -> IsRec(x) and IsBound(x.isScalarISA) and x.isScalarISA;

# register short variable names for SKLR(t) where t in  _SKLR_ISA_TYPES
# SKLR_i8, SKLR_ui8, ..., SKLR_f32, SKLR_f64, SKLR_f, SKLR_i
DoForAll( _SKLR_ISA_TYPES, t -> Assign(SKLR(t[1]).id(), SKLR(t[1])) );

# ISA.toScalar() returns scalar ISA of appropriate data type

SIMD_ISA.toScalar := self >> SKLR(self.gt());

ISA_Bridge.add(Class(CVT_SKLR_i16_i8,    ISA_Bridge_tcast, rec( props := [], isa_from := SKLR(T_Int(8)),  isa_to := SKLR(T_Int(16)) )));
ISA_Bridge.add(Class(CVT_SKLR_i16_ui8,   ISA_Bridge_tcast, rec( props := [], isa_from := SKLR(T_UInt(8)), isa_to := SKLR(T_Int(16)) )));
ISA_Bridge.add(Class(CVT_SKLR_ui16_ui8,  ISA_Bridge_tcast, rec( props := [], isa_from := SKLR(T_UInt(8)), isa_to := SKLR(T_UInt(16)) )));

ISA_Bridge.add(Class(CVT_SKLR_i32_i16,   ISA_Bridge_tcast, rec( props := [], isa_from := SKLR(T_Int(16)),  isa_to := SKLR(T_Int(32)) )));
ISA_Bridge.add(Class(CVT_SKLR_i32_ui16,  ISA_Bridge_tcast, rec( props := [], isa_from := SKLR(T_UInt(16)), isa_to := SKLR(T_Int(32)) )));
ISA_Bridge.add(Class(CVT_SKLR_ui32_ui16, ISA_Bridge_tcast, rec( props := [], isa_from := SKLR(T_UInt(16)), isa_to := SKLR(T_UInt(32)) )));

ISA_Bridge.add(Class(CVT_SKLR_i64_i32,   ISA_Bridge_tcast, rec( props := [], isa_from := SKLR(T_Int(32)),  isa_to := SKLR(T_Int(64)) )));
ISA_Bridge.add(Class(CVT_SKLR_i64_ui32,  ISA_Bridge_tcast, rec( props := [], isa_from := SKLR(T_UInt(32)), isa_to := SKLR(T_Int(64)) )));
ISA_Bridge.add(Class(CVT_SKLR_ui64_ui32, ISA_Bridge_tcast, rec( props := [], isa_from := SKLR(T_UInt(32)), isa_to := SKLR(T_UInt(64)) )));

ISA_Bridge.add(Class(CVT_SKLR_f32_i32,   ISA_Bridge_tcast, rec( props := [], isa_from := SKLR(T_Int(32)),  isa_to := SKLR(T_Real(32)) )));
ISA_Bridge.add(Class(CVT_SKLR_f32_ui32,  ISA_Bridge_tcast, rec( props := [], isa_from := SKLR(T_UInt(32)), isa_to := SKLR(T_Real(32)) )));

ISA_Bridge.add(Class(CVT_SKLR_f32_i64,   ISA_Bridge_tcast, rec( props := [], isa_from := SKLR(T_Int(64)),  isa_to := SKLR(T_Real(32)) )));
ISA_Bridge.add(Class(CVT_SKLR_f32_ui64,  ISA_Bridge_tcast, rec( props := [], isa_from := SKLR(T_UInt(64)), isa_to := SKLR(T_Real(32)) )));

ISA_Bridge.add(Class(CVT_SKLR_f64_i32,   ISA_Bridge_tcast, rec( props := [], isa_from := SKLR(T_Int(32)),  isa_to := SKLR(T_Real(64)) )));
ISA_Bridge.add(Class(CVT_SKLR_f64_ui32,  ISA_Bridge_tcast, rec( props := [], isa_from := SKLR(T_UInt(32)), isa_to := SKLR(T_Real(64)) )));

ISA_Bridge.add(Class(CVT_SKLR_f64_i64,   ISA_Bridge_tcast, rec( props := [], isa_from := SKLR(T_Int(64)),  isa_to := SKLR(T_Real(64)) )));
ISA_Bridge.add(Class(CVT_SKLR_f64_ui64,  ISA_Bridge_tcast, rec( props := [], isa_from := SKLR(T_UInt(64)), isa_to := SKLR(T_Real(64)) )));

ISA_Bridge.add(Class(CVT_SKLR_f64_f32,   ISA_Bridge_tcast, rec( props := [], isa_from := SKLR(T_Real(32)), isa_to := SKLR(T_Real(64)) )));

# helper class: integer to larger unsigned integer with saturation on zero
Class(ISA_Bridge_IxIntToUInt, ISA_Bridge_I, rec(
    props    := ["saturation"],
    code     := (self, y, x, opts) >> assign( self._y(y,0), tcast(self.isa_to.t, max(self.isa_from.t.zero(), self._x(x,0))) )
));

ISA_Bridge.add(Class(CVT_SKLR_ui16_i8_wrap,  ISA_Bridge_wrap,        rec( isa_from := SKLR(T_Int(8)),  isa_to := SKLR(T_UInt(16)))));
ISA_Bridge.add(Class(CVT_SKLR_ui16_i8_sat,   ISA_Bridge_IxIntToUInt, rec( isa_from := SKLR(T_Int(8)),  isa_to := SKLR(T_UInt(16)))));
ISA_Bridge.add(Class(CVT_SKLR_ui32_i16_wrap, ISA_Bridge_wrap,        rec( isa_from := SKLR(T_Int(16)), isa_to := SKLR(T_UInt(32)))));
ISA_Bridge.add(Class(CVT_SKLR_ui32_i16_sat,  ISA_Bridge_IxIntToUInt, rec( isa_from := SKLR(T_Int(16)), isa_to := SKLR(T_UInt(32)))));
ISA_Bridge.add(Class(CVT_SKLR_ui64_i32_wrap, ISA_Bridge_wrap,        rec( isa_from := SKLR(T_Int(32)), isa_to := SKLR(T_UInt(64)))));
ISA_Bridge.add(Class(CVT_SKLR_ui64_i32_sat,  ISA_Bridge_IxIntToUInt, rec( isa_from := SKLR(T_Int(32)), isa_to := SKLR(T_UInt(64)))));

# now from larger to smaller
# helper class: larger data type to smaller with max-min saturation
Class(ISA_Bridge_IxSat, ISA_Bridge_I, rec(
    props    := ["saturation"],
    code     := (self, y, x, opts) >> let(
         r := self.isa_to.t.range(),
         a := self.isa_from.t.value(r.min),
         b := self.isa_from.t.value(r.max),
         assign( self._y(y,0), tcast(self.isa_to.t, max(a, min(b, self._x(x,0))))) )
));

#NOTE: add SReduce for max(Value, a) where Value is leq a.t.range().min
#NOTE: add SReduce for max(Value, max(Value, a))
ISA_Bridge.add(Class(CVT_SKLR_i32_i64_wrap,   ISA_Bridge_wrap,  rec( isa_from := SKLR(T_Int(64)),  isa_to := SKLR(T_Int(32)) )));
ISA_Bridge.add(Class(CVT_SKLR_i32_i64_sat,    ISA_Bridge_IxSat, rec( isa_from := SKLR(T_Int(64)),  isa_to := SKLR(T_Int(32)) )));
ISA_Bridge.add(Class(CVT_SKLR_ui32_i64_wrap,  ISA_Bridge_wrap,  rec( isa_from := SKLR(T_Int(64)),  isa_to := SKLR(T_UInt(32)) )));
ISA_Bridge.add(Class(CVT_SKLR_ui32_i64_sat,   ISA_Bridge_IxSat, rec( isa_from := SKLR(T_Int(64)),  isa_to := SKLR(T_UInt(32)) )));
ISA_Bridge.add(Class(CVT_SKLR_i32_ui64_wrap,  ISA_Bridge_wrap,  rec( isa_from := SKLR(T_UInt(64)), isa_to := SKLR(T_Int(32)) )));
ISA_Bridge.add(Class(CVT_SKLR_i32_ui64_sat,   ISA_Bridge_IxSat, rec( isa_from := SKLR(T_UInt(64)), isa_to := SKLR(T_Int(32)) )));
ISA_Bridge.add(Class(CVT_SKLR_ui32_ui64_wrap, ISA_Bridge_wrap,  rec( isa_from := SKLR(T_UInt(64)), isa_to := SKLR(T_UInt(32)) )));
ISA_Bridge.add(Class(CVT_SKLR_ui32_ui64_sat,  ISA_Bridge_IxSat, rec( isa_from := SKLR(T_UInt(64)), isa_to := SKLR(T_UInt(32)) )));

ISA_Bridge.add(Class(CVT_SKLR_i16_i32_wrap,   ISA_Bridge_wrap,  rec( isa_from := SKLR(T_Int(32)),  isa_to := SKLR(T_Int(16)) )));
ISA_Bridge.add(Class(CVT_SKLR_i16_i32_sat,    ISA_Bridge_IxSat, rec( isa_from := SKLR(T_Int(32)),  isa_to := SKLR(T_Int(16)) )));
ISA_Bridge.add(Class(CVT_SKLR_ui16_i32_wrap,  ISA_Bridge_wrap,  rec( isa_from := SKLR(T_Int(32)),  isa_to := SKLR(T_UInt(16)) )));
ISA_Bridge.add(Class(CVT_SKLR_ui16_i32_sat,   ISA_Bridge_IxSat, rec( isa_from := SKLR(T_Int(32)),  isa_to := SKLR(T_UInt(16)) )));
ISA_Bridge.add(Class(CVT_SKLR_i16_ui32_wrap,  ISA_Bridge_wrap,  rec( isa_from := SKLR(T_UInt(32)), isa_to := SKLR(T_Int(16)) )));
ISA_Bridge.add(Class(CVT_SKLR_i16_ui32_sat,   ISA_Bridge_IxSat, rec( isa_from := SKLR(T_UInt(32)), isa_to := SKLR(T_Int(16)) )));
ISA_Bridge.add(Class(CVT_SKLR_ui16_ui32_wrap, ISA_Bridge_wrap,  rec( isa_from := SKLR(T_UInt(32)), isa_to := SKLR(T_UInt(16)) )));
ISA_Bridge.add(Class(CVT_SKLR_ui16_ui32_sat,  ISA_Bridge_IxSat, rec( isa_from := SKLR(T_UInt(32)), isa_to := SKLR(T_UInt(16)) )));

ISA_Bridge.add(Class(CVT_SKLR_i8_i16_wrap,    ISA_Bridge_wrap,  rec( isa_from := SKLR(T_Int(16)),  isa_to := SKLR(T_Int(8)) )));
ISA_Bridge.add(Class(CVT_SKLR_i8_i16_sat,     ISA_Bridge_IxSat, rec( isa_from := SKLR(T_Int(16)),  isa_to := SKLR(T_Int(8)) )));
ISA_Bridge.add(Class(CVT_SKLR_ui8_i16_wrap,   ISA_Bridge_wrap,  rec( isa_from := SKLR(T_Int(16)),  isa_to := SKLR(T_UInt(8)) )));
ISA_Bridge.add(Class(CVT_SKLR_ui8_i16_sat,    ISA_Bridge_IxSat, rec( isa_from := SKLR(T_Int(16)),  isa_to := SKLR(T_UInt(8)) )));
ISA_Bridge.add(Class(CVT_SKLR_i8_ui16_wrap,   ISA_Bridge_wrap,  rec( isa_from := SKLR(T_UInt(16)), isa_to := SKLR(T_Int(8)) )));
ISA_Bridge.add(Class(CVT_SKLR_i8_ui16_sat,    ISA_Bridge_IxSat, rec( isa_from := SKLR(T_UInt(16)), isa_to := SKLR(T_Int(8)) )));
ISA_Bridge.add(Class(CVT_SKLR_ui8_ui16_wrap,  ISA_Bridge_wrap,  rec( isa_from := SKLR(T_UInt(16)), isa_to := SKLR(T_UInt(8)) )));
ISA_Bridge.add(Class(CVT_SKLR_ui8_ui16_sat,   ISA_Bridge_IxSat, rec( isa_from := SKLR(T_UInt(16)), isa_to := SKLR(T_UInt(8)) )));

ISA_Bridge.add(Class(CVT_SKLR_f32_f64_trunc,   ISA_Bridge_wrap, rec( props := [], isa_from := SKLR(T_Real(64)), isa_to := SKLR(T_Real(32)) )));

Class(ISA_Bridge_trunc_wrap, ISA_Bridge_I, rec(
    range := (self) >> self.isa_to.t.range(),
    props := ["trunc", "wraparound"],
    code  := (self, y, x, opts) >> assign( self._y(y,0), tcast(self.isa_to.t, self._x(x,0)) ),
));

ISA_Bridge.add(Class(CVT_SKLR_i32_f32_trunc,   ISA_Bridge_trunc_wrap, rec( isa_from := SKLR(T_Real(32)), isa_to := SKLR(T_Int(32))  )));
ISA_Bridge.add(Class(CVT_SKLR_ui32_f32_trunc,  ISA_Bridge_trunc_wrap, rec( isa_from := SKLR(T_Real(32)), isa_to := SKLR(T_UInt(32)) )));

ISA_Bridge.add(Class(CVT_SKLR_i32_f64_trunc,   ISA_Bridge_trunc_wrap, rec( isa_from := SKLR(T_Real(32)), isa_to := SKLR(T_Int(64))  )));
ISA_Bridge.add(Class(CVT_SKLR_ui32_f64_trunc,  ISA_Bridge_trunc_wrap, rec( isa_from := SKLR(T_Real(32)), isa_to := SKLR(T_UInt(64)) )));

ISA_Bridge.add(Class(CVT_SKLR_i64_f32_trunc,   ISA_Bridge_trunc_wrap, rec( isa_from := SKLR(T_Real(64)), isa_to := SKLR(T_Int(32))  )));
ISA_Bridge.add(Class(CVT_SKLR_ui64_f32_trunc,  ISA_Bridge_trunc_wrap, rec( isa_from := SKLR(T_Real(64)), isa_to := SKLR(T_UInt(32)) )));

ISA_Bridge.add(Class(CVT_SKLR_i64_f64_trunc,   ISA_Bridge_trunc_wrap, rec( isa_from := SKLR(T_Real(64)), isa_to := SKLR(T_Int(64))  )));
ISA_Bridge.add(Class(CVT_SKLR_ui64_f64_trunc,  ISA_Bridge_trunc_wrap, rec( isa_from := SKLR(T_Real(64)), isa_to := SKLR(T_UInt(64)) )));

Class(ISA_Bridge_round_wrap, ISA_Bridge_I, rec(
    range := (self) >> self.isa_to.t.range(),
    props := ["round", "wraparound"],
    code  := (self, y, x, opts) >> assign( self._y(y,0), tcast(self.isa_to.t, add(self._x(x,0), self.isa_from.t.value(0.5))) ),
));

ISA_Bridge.add(Class(CVT_SKLR_i32_f32_round,   ISA_Bridge_round_wrap, rec( isa_from := SKLR(T_Real(32)), isa_to := SKLR(T_Int(32))  )));
ISA_Bridge.add(Class(CVT_SKLR_ui32_f32_round,  ISA_Bridge_round_wrap, rec( isa_from := SKLR(T_Real(32)), isa_to := SKLR(T_UInt(32)) )));

ISA_Bridge.add(Class(CVT_SKLR_i32_f64_round,   ISA_Bridge_round_wrap, rec( isa_from := SKLR(T_Real(32)), isa_to := SKLR(T_Int(64))  )));
ISA_Bridge.add(Class(CVT_SKLR_ui32_f64_round,  ISA_Bridge_round_wrap, rec( isa_from := SKLR(T_Real(32)), isa_to := SKLR(T_UInt(64)) )));

ISA_Bridge.add(Class(CVT_SKLR_i64_f32_round,   ISA_Bridge_round_wrap, rec( isa_from := SKLR(T_Real(64)), isa_to := SKLR(T_Int(32))  )));
ISA_Bridge.add(Class(CVT_SKLR_ui64_f32_round,  ISA_Bridge_round_wrap, rec( isa_from := SKLR(T_Real(64)), isa_to := SKLR(T_UInt(32)) )));

ISA_Bridge.add(Class(CVT_SKLR_i64_f64_round,   ISA_Bridge_round_wrap, rec( isa_from := SKLR(T_Real(64)), isa_to := SKLR(T_Int(64))  )));
ISA_Bridge.add(Class(CVT_SKLR_ui64_f64_round,  ISA_Bridge_round_wrap, rec( isa_from := SKLR(T_Real(64)), isa_to := SKLR(T_UInt(64)) )));

ISA_Bridge.add(Class(CVT_SKLR_f32_clip32i, ISA_Bridge_I, rec(
    isa_from    := SKLR(T_Real(32)),
    isa_to      := SKLR(T_Real(32)),
    props       := ["saturation"],
    range       := (self) >> RangeT(self.clip.min, self.clip.max, T_Real(32).range().eps),
    clip        := rec( min := T_Int(32).range().min, max := T_Int(32).range().max ),
    code := (self, y, x, opts) >> assign( self._y(y,0), min(max(self._x(x,0), self.isa_from.t.value(self.clip.min)), self.isa_from.t.value(self.clip.max)) ),
)));

ISA_Bridge.add(Class(CVT_SKLR_f32_clip32ui, CVT_SKLR_f32_clip32i, rec(
    clip        := rec( min := T_UInt(32).range().min, max := T_UInt(32).range().max ),
)));

ISA_Bridge.add(Class(CVT_SKLR_f32_clip16i, CVT_SKLR_f32_clip32i, rec(
    clip        := rec( min := T_Int(16).range().min,  max := T_Int(16).range().max ),
)));

ISA_Bridge.add(Class(CVT_SKLR_f32_clip16ui, CVT_SKLR_f32_clip32i, rec(
    clip        := rec( min := T_UInt(16).range().min, max := T_UInt(16).range().max ),
)));

ISA_Bridge.add(Class(CVT_SKLR_f32_clip8i, CVT_SKLR_f32_clip32i, rec(
    clip        := rec( min := T_Int(8).range().min,  max := T_Int(8).range().max ),
)));

ISA_Bridge.add(Class(CVT_SKLR_f32_clip8ui, CVT_SKLR_f32_clip32i, rec(
    clip        := rec( min := T_UInt(8).range().min, max := T_UInt(8).range().max ),
)));


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

Declare(TDST2, TDST3, TDST4);

_DST_CONST := 2.5;

#######################################################################################
#   tSPL DST rules

#F DST4(<n>) - Discrete Sine Transform, Type IV, non-terminal
#F Definition: (n x n)-matrix [ sin((k-1/2)*(l-1/2)*pi/n) | k,l = 1...n ]
#F Note:       DST4 is symmetric
#F Example:    DST4(8)
Class(TDST4, TaggedNonTerminal, rec(
    abbrevs := [ N -> Checked(IsInt(N), N >= 1, [N]) ] ,
    dims := self >> [self.params[1], self.params[1]],
    isReal := True,
    terminate := self >> Mat(DST_IVunscaled(self.params[1])),
    transpose := self >> Copy(self),
    SmallRandom := () -> Random([2, 3, 4, 6, 8, 9, 12, 16, 18, 24, 27, 32]),
    normalizedArithCost := (self) >> let(n := self.params[1], IntDouble(_DST_CONST * n * d_log(n) / d_log(2)))
));

NewRulesFor(TDST4, rec(
    DST4_CT_tSPL := rec(
        switch := false,

        applicable := (self, t) >> let(
            P:=t.params,
            P[1] > 2
            and ForAny(DivisorPairs(2*P[1]), d ->
                IsEvenInt(d[1]) and IsEvenInt(d[2])
            )
            and t.hasTags()
        ),

        children := (self, t) >> let(
            tags := t.getTags(),
            N := t.params[1],
            Map2(
                Filtered(DivisorPairs(2*N), d -> IsEvenInt(d[1]) and IsEvenInt(d[2])),
                (m,n) -> List([
                    TPrm(IJ(N, n/2)),
                    TTensorI(condM(m,m/2) * PRDFT3(m,-1).transpose(), n/2, AVec, AVec),
                    TTensorI(L(n, 2) * PRDFT3(n) * condK(n, 2), m/2, APar, AVec),
                    TPrm(IJ(N, m/2))
                ], i -> i.setTags(tags))
            )
        ),

        apply := (self, t, C, Nonterms) >> let(
            N:=t.params[1],
            n:=2*Rows(Nonterms[1].params[1].params[2]),
            m:=2*Rows(Nonterms[4].params[1].params[2]),
            Grp(
                C[1] * C[2]
                * ConjDiag(RC(Diag(fPrecompute(diagMul(fConst(TComplex, N/2, -E(4)),
                        fCompose(dOmega(8 * N, 1),
                            diagTensor(dLin(N/m, 2, 1, TInt), dLin(m/2, 2, 1, TInt))))
                ))), L(N, m), L(N, n/2))
            )
            * C[3] * C[4]
        )
    )
));

#F TDST2(<n>) - Discrete Sine Transform, Type II, non-terminal
#F Definition: (n x n)-matrix [ sin(k*(l+1/2)*pi/n) | k,l = 0...n-1 ]
#F Note:       DST2 is the transpose of DST3
#F Example:    DST2(8)
Class(TDST2, TaggedNonTerminal, rec(
    abbrevs := [ N -> Checked(IsInt(N), N >= 1, [N]) ] ,
    dims := self >> [self.params[1], self.params[1]],
    isReal := True,
    terminate := self >> Mat(DST_IIunscaled(self.params[1])),
    transpose := self >> TDST3(self.params[1]),
    SmallRandom := () -> Random([2,3,4,5,6,8,9,10,12,15,16,18,24,27,30,32]),
    normalizedArithCost := (self) >> let(n := self.params[1], IntDouble(_DST_CONST * n * d_log(n) / d_log(2)))
));

NewRulesFor(TDST2, rec(
    DST2_DST4_tSPL := rec(
        switch := false,
        applicable := (self, t) >> true,
        children := (self, t) >> [[ TS(t.params[1]).withTags(t.getTags()), TDST4(t.params[1]).withTags(t.getTags()) ]],
        apply := (self, t, C, Nonterms) >> let(P := t.params[1],
            Diag(Concat(List([1..P-1], i->V(1.0)), [V(2.0)])) *
            C[1] *
            C[2] *
            Diag(List([0..P - 1], i -> 1/(2 * CosPi((2*i + 1)/(4*P)))))
        )
    ))
);

#F TDST3(<n>) - Discrete Sine Transform, Type III, non-terminal
#F Definition: (n x n)-matrix [ sin((k-1/2)*l*pi/n) | k,l = 1...n ]
#F Note:       DST3 is the transpose of DST2
#F Example:    DST3(8)
#F Scaled variant (not supported) is:
#F                [ a_l*sin(k*(l-1/2)*pi/n) | k,l = 1...n ]
#F            with  a_j = 1/sqrt(2) for j = n and = 1 else.
Class(TDST3, TaggedNonTerminal, rec(
    abbrevs := [ N -> Checked(IsInt(N), N >= 1, [N]) ] ,
    dims := self >> [self.params[1], self.params[1]],
    isReal := True,
    terminate := self >> Mat(DST_IIIunscaled(self.params[1])),
    transpose := self >> TDST2(self.params[1]),
    SmallRandom := () -> Random([2,3,4,5,6,8,9,10,12,15,16,18,24,27,30,32]),
    normalizedArithCost := (self) >> let(n := self.params[1], IntDouble(_DST_CONST * n * d_log(n) / d_log(2)))
));

NewRulesFor(TDST3, rec(
    DST3_DST4_tSPL := rec(
        switch := false,
        applicable := (self, t) >> true,
        children := (self, t) >> [[ TDST4(t.params[1]).withTags(t.getTags()), TS(t.params[1]).withTags(t.getTags()) ]],
        apply := (self, t, C, Nonterms) >> let(P := t.params[1],
            Diag(List([0..P - 1], i -> 1/(2 * CosPi((2*i + 1)/(4*P))))) *
            C[1] * C[2].transpose() *
            Diag(Concat(List([1..P-1], i->V(1.0)), [V(2.0)])))
    )
));


text:
# You can't break an outer loop unless you return from the whole function.
n := 40;
a := List([1 .. n], i -> List([1 .. n], j -> Random(1, 20)));;

Find := function(a, x)
    local i, j, n;
    n := Length(a);
    for i in [1 .. n] do
        for j in [1 .. n] do
            if a[i][j] = x then
                return [i, j];
            fi;
        od;
    od;
    return fail;
end;

Find(a, 20);

max_stars_repo_name: mullikine/RosettaCodeData