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

#####################
# General Rule
#####################

PRF34_CT_Children := (N,k,DFTp,PRF3,PRFt) -> Map2(DivisorPairs(N),
    (m,n) -> When(IsEvenInt(n),
	[ DFTp(m, k), PRF3(n, k) ],
	[ DFTp(m, k), PRF3(n, k), PRFt(m, k) ])
);

PRF34_CT_Rule := (N,k,C,Conj,Tw) -> let(
    m:=Rows(C[1]), n:=Cols(C[2]), Nc:=Int((N+1)/2),
    nf:=Int(n/2), nc:=Int((n+1)/2), mf:=Int(m/2), mc:=Int((m+1)/2), j:=Ind(nf),

    SUM(
	ISum(j, RC(Scat(BH(Nc, N-1, m, j, n))) *
	        Conj * RC(C[1]) * Tw(j) *
		RC(Gath(H(nc*m, m, j, nc)))
	),
	When(IsEvenInt(n), [],
	     RC(Scat(H(Nc, mc, nf, n))) * C[3] * Gath(H(2*nc*m, m, 2*nf, 2*nc)))
    ) *
    Tensor(I(m), C[2]) *
    L(N,m)
);

IPRF34_CT_Rule := (N,k,C,Conj,Tw) -> let(
    m:=Rows(C[1]), n:=Rows(C[2]), Nc:=Int((N+1)/2),
    nf:=Int(n/2), nc:=Int((n+1)/2), mf:=Int(m/2), mc:=Int((m+1)/2), j:=Ind(nf),

    L(N,n) *
    Tensor(I(m), C[2]) *
    SUM(
	ISum(j,
	    RC(Scat(H(nc*m, m, j, nc))) *
	    Tw(j) *
	    RC(C[1]) *
	    Conj *
	    RC(Gath(BH(Nc, N-1, m, j, n)))
	),
	When(IsEvenInt(n), [],
	     Scat(H(2*nc*m, m, 2*nf, 2*nc)) * C[3] * RC(Gath(H(Nc, mc, nf, n))))
    )
);

#####################
# Special Cases
#####################

RulesFor(PRDFT4, rec(
    PRDFT4_Base1 := BaseRule(PRDFT4, [1, @]),
    PRDFT4_Base2 := rec(
	isApplicable := P -> P[1]=2,
	rule := (P, C) -> F(2) * Mat(1/2*MatSPL(PDHT4(2,P[2])))  # DHT4(2) is a diagonal
    ),
    PRDFT4_CT := rec(
	isApplicable := P -> not IsPrime(P[1]),
	allChildren  := P -> PRF34_CT_Children(P[1], P[2], DFT2, PRDFT3, PRDFT4),
	rule := (P,C) -> let(N:=P[1], k:=P[2], m:=Rows(C[1]),
	    PRF34_CT_Rule(N, k, C, Diag(BHD(m,-1,1)), j->RC(Diag(fPrecompute(Twid(N,m,k,1/2,1/2,j))))))),
    PRDFT4_Trig := rec(
	isApplicable := P -> P[1] > 2 and IsEvenInt(P[1]) and P[2]=1,
	allChildren := P -> [[ DCT4(P[1]/2), DST4(P[1]/2) ]],
	rule := (P, C) -> let(n:=P[1]/2,
	    L(2*n, n) * DirectSum(C[1], C[2]) * SymSplit4(n)))
));

RulesFor(IPRDFT4, rec(
    IPRDFT4_Base1 := BaseRule(IPRDFT4, [1, @]),
    IPRDFT4_Base2 := rec(
	isApplicable := P -> P[1]=2,
	rule := (P, C) -> Mat(1/2*MatSPL(PDHT4(2,P[2]))) * F(2)  # DHT4(2) is a diagonal
    ),
    IPRDFT4_CT := rec(
	isApplicable := P -> not IsPrime(P[1]),
	allChildren  := P -> PRF34_CT_Children(P[1], P[2], DFT3, IPRDFT2, IPRDFT4),
	rule := (P,C) -> let(N:=P[1], k:=P[2], m:=Rows(C[1]),
	    IPRF34_CT_Rule(N, k, C, Diag(BHD(m,-1,1)), j->RC(Diag(fPrecompute(Twid(N,m,k,1/2,1/2,j)))))))
));

RulesFor(PRDFT3, rec(
    PRDFT3_Base1 := BaseRule(PRDFT3, [1, @]),
    PRDFT3_Base2 := rec(
	isApplicable := P -> P[1]=2,
	rule := (P, C) -> Cond(
	    P[2] mod 4 = 1, I(2),
	    P[2] mod 4 = 3, Diag(1,-1),
	    Error("Bad second parameter for PRDFT3(n,k), gcd(n,k)=1 does not hold"))
    ),
    PRDFT3_CT := rec(
	isApplicable := P -> not IsPrime(P[1]),
	allChildren  := P -> PRF34_CT_Children(P[1], P[2], DFT1, PRDFT3, PRDFT3),
	rule := (P,C) -> let(N:=P[1], k:=P[2], m:=Rows(C[1]),
	    PRF34_CT_Rule(N, k, C, Diag(BHD(m,1,-1)), j->RC(Diag(fPrecompute(Twid(N,m,k,1/2,0,j))))))),
    PRDFT3_Trig := rec(
	isApplicable := P -> P[1] > 2 and IsEvenInt(P[1]) and P[2]=1,
	allChildren := P -> [[ DCT3(P[1]/2), DST3(P[1]/2) ]],
	rule := (P, C) -> let(n:=P[1]/2,
	    L(2*n, n) * DirectSum(C[1], C[2]) * SymSplit3(n))),

    PRDFT3_OddToPRDFT1 := rec(
	isApplicable := P -> IsOddInt(P[1]),
	allChildren := P -> [[ PRDFT(P[1], P[2]) ]],
	rule := (P, C) -> rperm_ev(P[1]) * C[1] * pdiag(P[1], P[2])
    )
));

RulesFor(IPRDFT2, rec(
    IPRDFT2_Base1 := BaseRule(IPRDFT2, [1, @]),
    IPRDFT2_Base2 := rec(
	isApplicable := P -> P[1]=2,
	rule := (P, C) -> Cond(
	    P[2] mod 4 = 1, Diag(2,-2),
	    P[2] mod 4 = 3, Diag(2,2),
	    Error("Bad second parameter for IPRDFT2(n,k), gcd(n,k)=1 does not hold"))
    ),
    IPRDFT2_CT := rec(
	isApplicable := P -> not IsPrime(P[1]),
	allChildren  := P -> PRF34_CT_Children(P[1], P[2], DFT1, IPRDFT2, IPRDFT2),
	rule := (P,C) -> let(N:=P[1], k:=P[2], m:=Rows(C[1]),
	    IPRF34_CT_Rule(N, k, C, Diag(BHD(m,1,-1)), j->RC(Diag(fPrecompute(Twid(N,m,k,1/2,0,j)))))))
));

RulesFor(PDHT4, rec(
    PDHT4_Base2 := BaseRule(PDHT4, [2, @]),
    PDHT4_CT := rec(
	isApplicable := P -> not IsPrime(P[1]),
	allChildren  := P -> PRF34_CT_Children(P[1], P[2], DFT2, PDHT3, PDHT4),
	rule := (P,C) -> let(N:=P[1], k:=P[2], m:=Rows(C[1]),
	    PRF34_CT_Rule(N, k, C, TopHalf(m, -J(2)),
		          j -> RCDiag(RCData(fPrecompute(Twid(N,m,-k,1/2,1/2,j))), -J(2))))),
    PDHT4_Trig := rec(
	isApplicable := P -> P[1] > 2 and IsEvenInt(P[1]) and P[2]=1,
	allChildren := P -> [[ DCT4(P[1]/2), DST4(P[1]/2) ]],
	rule := (P, C) -> let(n:=P[1]/2,
	    Tensor(I(n), F(2)) * L(2*n,n) * DirectSum(C[1], C[2]) * SymSplit4(n)))
));

RulesFor(PDHT3, rec(
    PDHT3_Base2 := rec(
	isApplicable := P -> P[1]=2,
	rule := (P, C) -> Cond(
	    P[2] mod 4 = 1, F(2),
	    P[2] mod 4 = 3, J(2)*F(2),
	    Error("Bad second parameter for PRDFT3(n,k), gcd(n,k)=1 does not hold"))
    ),

    PDHT3_CT := rec(
	isApplicable := P -> not IsPrime(P[1]),
	allChildren  := P -> PRF34_CT_Children(P[1], P[2], DFT1, PDHT3, PDHT3),
	rule := (P,C) -> let(N:=P[1], k:=P[2], m:=Rows(C[1]),
	    PRF34_CT_Rule(N, k, C, TopHalf(m, J(2)),
		          j -> RCDiag(RCData(fPrecompute(Twid(N,m,-k,1/2,0,j))), J(2))))),
    PDHT3_Trig := rec(
	isApplicable := P -> P[1] > 2 and IsEvenInt(P[1]) and P[2]=1,
	allChildren := P -> [[ DCT3(P[1]/2), DST3(P[1]/2) ]],
	rule := (P, C) -> let(n:=P[1]/2,
	    Tensor(I(n), F(2)) * L(2*n,n) * DirectSum(C[1], C[2]) * SymSplit3(n)))

));


text:
Anagrams := function(name)
  local f, p, L, line, word, words, swords, res, cur, r;
  words := [ ];
  swords := [ ];
  f := InputTextFile(name);
  while true do
    line := ReadLine(f);
    if line = fail then
      break;
    else
      word := Chomp(line);
      Add(words, word);
      Add(swords, SortedList(word));
    fi;
  od;
  CloseStream(f);
  p := SortingPerm(swords);
  L := Permuted(words, p);
  r := "";
  cur := [ ];
  res := [ ];
  for word in L do
    if SortedList(word) = r then
      Add(cur, word);
    else
      if Length(cur) > 0 then
        Add(res, cur);
      fi;
      r := SortedList(word);
      cur := [ word ];
    fi;
  od;
  if Length(cur) > 0 then
    Add(res, cur);
  fi;
  return Filtered(res, v -> Length(v) > 1);
end;

ana := Anagrams("my/gap/unixdict.txt");;

# What is the longest anagram sequence ?
Maximum(List(ana, Length));
# 5

# Which are they ?
Filtered(ana, v -> Length(v) = 5);
# [ [ "abel", "able", "bale", "bela", "elba" ],
#   [ "caret", "carte", "cater", "crate", "trace" ],
#   [ "angel", "angle", "galen", "glean", "lange" ],
#   [ "alger", "glare", "lager", "large", "regal" ],
#   [ "elan", "lane", "lean", "lena", "neal" ],
#   [ "evil", "levi", "live", "veil", "vile" ] ]

ext: gap

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


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

_dbgPrintRuleName := rule -> Cond(
    IsBound(rule.owner),
        Print(spiral._color(Blue, rule.name), "   /  ",
	      spiral._color(DarkBlue, rule.owner), "\n"),
    Print(spiral._color(Blue, rule.name), "\n")
);

_dbgRuleBelongsTo := (rule, rset) ->
    IsBound(rule.owner) and (rule.owner in rset);

#F DebugRewriting(<rsets.)
DebugRewriting := function(rsets)
    if rsets=true then
        rewrite.RuleTrace := rule -> _dbgPrintRuleName(rule);
        rewrite.RuleStatus := (rule, hd, str) -> Print(
	    spiral._color(Yellow, hd), ApplyFunc(Print, str));
    elif rsets=false then
        rewrite.RuleTrace := Ignore;
        rewrite.RuleStatus := Ignore;
    else
		if not IsList(rsets) then
			rsets := [rsets];
		fi;
        rewrite.RuleTrace := rule ->
			When(_dbgRuleBelongsTo(rule, rsets), _dbgPrintRuleName(rule));
        rewrite.RuleStatus := (rule, hd, str) ->
			When(_dbgRuleBelongsTo(rule, rsets),
		Print(spiral._color(Yellow, hd), ApplyFunc(Print, str)));
    fi;
end;

DebugRuleStrategies := function(switch)
    Constraint(IsBool(switch));
    if switch then
        rewrite.RuleTrace := rule -> _dbgPrintRuleName(rule);
        rewrite.RuleStrategyTrace := (i, rset, expr) -> Print(
			spiral._color(Red, i),
			spiral._color(Red, " ----------------\n"),
			rset, "\n",
			Doc(rset),
			expr, "\n");
    else
        rewrite.RuleTrace := Ignore;
        rewrite.RuleStrategyTrace := Ignore;
    fi;
end;

ext: gi


text:
# The Grassmann graph J_q(n, d) of d-dimensional subspaces of F_q^n.
BindGlobal("GrassmannGraph", function(q, n, d)
    local G;
    G := SubspaceGraph(GL(n, q), Elements, GF(q)^n, d, true);
    G.duality := GrassmannDualityFunction;
    G.primality := GrassmannDualityFunction;
    return G;
end);

# The doubled Grassmann graph 2J_q(2d+1, d) of d- and (d+1)-dimensional
# subspaces of F_q^{2d+1}.
BindGlobal("DoubledGrassmannGraph", function(q, d)
    local n, G, V, dp;
    n := 2*d+1;
    V := GF(q)^n;
    dp := DirectProduct(GL(n, q), Group((1, 2)));
    G := Graph(dp, Union(Subspaces(V, d), Subspaces(V, d+1)),
                OnDoubledGrassmann(V, dp), SymmetrizedInclusion, true);
    G.halfDuality := GrassmannDualityFunction;
    G.halfPrimality := GrassmannDualityFunction;
    return G;
end);

# The twisted Grassmann graph TG_d(q) of (d+1)-dimensional subspaces of
# F_q^{2d+1} which are not subspaces of a hyperplane H, and (d-1)-dimensional
# subspaces of H.
BindGlobal("TwistedGrassmannGraph", function(d, q)
    local F, H, S, V, n;
    n := 2*d+1;
    F := GF(q)^n;
    H := Subspace(F, CanonicalBasis(F){[1..n-1]}, "basis");
    V := Union(Filtered(Subspaces(F, d+1), x -> not IsSubset(H, x)),
            Subspaces(H, d-1));
    S := OnSubspaces(F);
    return Graph(Stabilizer(GL(n, q), H, S), V, S, function(x, y)
        local dx;
        dx := Dimension(x);
        if dx = Dimension(y) then
            return Dimension(Intersection(x, y)) = dx-1;
        else
            return IsSubset(x, y) or IsSubset(y, x);
        fi;
    end, true);
end);

# The polar graph O^{(+/-)}(d, q) of isotropic lines of F_q^d
# with respect to a nondegenerate quadratic form.
BindGlobal("PolarGraphO", function(arg)
    local d, e, q, B, G, H, Q, V;
    if Length(arg) < 2 then
        Error("at least two arguments expected");
        return fail;
    elif Length(arg) = 2 then
        e := 0;
        d := arg[1];
        q := arg[2];
    else
        e := arg[1];
        d := arg[2];
        q := arg[3];
    fi;
    G := GO(e, d, q);
    Q := InvariantQuadraticForm(G).matrix;
    B := Q + TransposedMat(Q);
    V := GF(q)^d;
    H := Graph(G, IsotropicSpacesQuadraticForm(Q)(Subspaces(V, 1)),
                OnSubspaces(V), function(x, y)
                    return x <> y and IsZero(Elements(x)[2]*B*Elements(y)[2]);
                end, true);
    H.duality := Sum;
    H.primality := Intersection;
    return H;
end);

# The polar graph NO^{+/-}orth(d, q) of nonisotropic lines of F_q^d
# with respect to a nondegenerate quadratic form.
BindGlobal("PolarGraphNOorth", function(e, d, q)
    local f, z, B, G, H, Q, V;
    if d mod 2 = 1 then
        f := 0;
        if q mod 4 = 3 and d mod 4 = 1 then
            e := -e;
        fi;
    else
        f := e;
    fi;
    G := GO(f, d, q);
    Q := InvariantQuadraticForm(G).matrix;
    B := Q + TransposedMat(Q);
    V := GF(q)^d;
    z := Z(q)^((1-e)/2);
    H := Graph(G, NonisotropicSpacesQuadraticForm(Q, z)(Subspaces(V, 1)),
                OnSubspaces(V), function(x, y)
                    return x <> y and IsZero(Elements(x)[2]*B*Elements(y)[2]);
                end, true);
    return H;
end);

# The polar graph Sp(d, q) of isotropic lines of F_q^d
# with respect to a nondegenerate symplectic form.
BindGlobal("PolarGraphSp", function(d, q)
    local G, H, Q, V;
    G := Sp(d, q);
    Q := InvariantBilinearForm(G).matrix;
    V := GF(q)^d;
    H := Graph(G, IsotropicSpacesQuadraticForm(Q)(Subspaces(V, 1)),
                OnSubspaces(V), function(x, y)
                    return x <> y and IsZero(Elements(x)[2]*Q*Elements(y)[2]);
                end, true);
    H.duality := Sum;
    H.primality := Intersection;
    return H;
end);

# The polar graph U(d, r) of isotropic lines of F_{r^2}^d
# with respect to a nondegenerate Hermitean form.
BindGlobal("PolarGraphU", function(d, r)
    local c, B, F, G, H, Q, V;
    G := GU(d, r);
    Q := InvariantSesquilinearForm(G).matrix;
    V := GF(r^2)^d;
    B := Elements(CanonicalBasis(V));
    c := Conjugates(GF(r^2), GF(r), Z(r^2));
    F := x -> List(x, y -> y^r);
    H := Graph(G, IsotropicSpacesSesquilinearForm(Q, r)(Subspaces(V, 1)),
            OnSubspaces(V), function(x, y)
                return x <> y and IsZero(Elements(x)[2]*Q*F(Elements(y)[2]));
            end, true);
    H.duality := Sum;
    H.primality := Intersection;
    return H;
end);

# The dual polar graph B_d(q) of isotropic d-dimensional subspaces of
# F_q^{2d+1} with respect to a nondegenerate quadratic form.
BindGlobal("DualPolarGraphB", function(d, q)
    local F, G, S, e, z, invt;
    e := 2*d+1;
    F := GF(q)^e;
    G := GO(e, q);
    if q mod 2 = 1 then
        invt := true;
        S := IsotropicSpacesQuadraticForm(InvariantQuadraticForm(G).matrix);
    else
        invt := false;
        S := [Subspace(F, Elements(CanonicalBasis(F)){[2..d+1]}, "basis")];
    fi;
    return SubspaceGraph(G, S, F, d, invt);
end);

# The dual polar graph C_d(q) of isotropic d-dimensional subspaces of
# F_q^{2d} with respect to a nondegenerate symplectic form.
BindGlobal("DualPolarGraphC", function(d, q)
    local F, e;
    e := 2*d;
    F := GF(q)^e;
    return SubspaceGraph(Sp(e, q),
        [Subspace(F, Elements(CanonicalBasis(F)){[1..d]}, "basis")],
        F, d, false);
end);

# The dual polar graph D_d(q) of isotropic d-dimensional subspaces of
# F_q^{2d} with respect to a nondegenerate quadratic form of Witt index d.
BindGlobal("DualPolarGraphD", function(d, q)
    local F, G, S, e, invt;
    e := 2*d;
    F := GF(q)^e;
    G := GO(1, e, q);
    if q mod 2 = 1 then
        invt := true;
        S := IsotropicSpacesQuadraticForm(InvariantQuadraticForm(G).matrix);
    else
        invt := false;
        S := [Subspace(F, Elements(CanonicalBasis(F)){[1,3..e-1]}, "basis")];
    fi;
    return SubspaceGraph(G, S, F, d, invt);
end);

# The dual polar graph ^2D_{d+1}(q) of isotropic d-dimensional subspaces of
# F_q^{2d+2} with respect to a nondegenerate quadratic form of Witt index d.
BindGlobal("DualPolarGraph2D", function(d, q)
    local G, e;
    e := 2*d+2;
    G := GO(-1, e, q);
    return SubspaceGraph(G,
        IsotropicSpacesQuadraticForm(InvariantQuadraticForm(G).matrix),
        GF(q)^e, d);
end);

# The dual polar graph ^2A_{e-1}(r) of isotropic [e/2]-dimensional
# subspaces of F_{r^2}^e with respect to a nondegenerate Hermitean form.
BindGlobal("DualPolarGraph2A", function(e, r)
    local B, F, c, d;
    F := GF(r^2)^e;
    B := Elements(CanonicalBasis(F));
    c := Conjugates(GF(r^2), GF(r), Z(r^2));
    d := Int(e/2);
    return SubspaceGraph(GU(e, r),
        [Subspace(F, B{[1..d]} + (c[1]-c[2])*B{e+1-[1..d]}, "basis")],
        F, d, false);
end);

# The Doro graph of nonisotropic 1-dimensional subspaces of F_q^4 with respect
# to a nondegenerate quadratic form. It is distance-regular for q = 4, 5.
BindGlobal("DoroGraph", function(q)
    local G, V;
    V := GF(q)^4;
    G := GO(-1, 4, q);
    return Graph(G, [Subspace(V, BasisVectors(Basis(V)){[4]}, "basis")],
        OnSubspaces(V), IsHyperbolic(InvariantQuadraticForm(G).matrix));
end);

# The unitary nonisotropics graph of 1-dimensional subspaces of F_(r^2)^3 with
# respect to a nondegenerate sesquilinear form.
BindGlobal("UnitaryNonisotropicsGraph", function(r)
    local G, V;
    V := GF(r^2)^3;
    G := GU(3, r);
    return Graph(G, [Subspace(V, BasisVectors(Basis(V)){[2]}, "basis")],
        OnSubspaces(V), IsOrthogonal(InvariantSesquilinearForm(G).matrix, r));
end);

max_stars_repo_name: jaanos/gap-graphs

max_issues_repo_name: jaanos/gap-graphs

max_forks_repo_name: jaanos/gap-graphs


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

_neon_binop_ev := rec(
    ev := self >> self.t.value(self.semantic(self._vval(1), self._vval(2), [])).ev(),
);

############### HALF NEON #############
Class(vunpacklo_half, _neon_binop_ev, VecExp_2.binary(), rec(
    semantic := (in1, in2, p) -> unpacklo(in1, in2, 2, 1)));

Class(vunpackhi_half, _neon_binop_ev, VecExp_2.binary(), rec(
    semantic := (in1, in2, p) -> unpackhi(in1, in2, 2, 1)));

Class(vrev_half, VecExp_2.unary(), rec(
    semantic := (in1, p) -> vrev64(in1,2),
));

Class(vtrnq_half, _neon_binop_ev, VecExp_2.binary(), rec(
    semantic := (in1, in2, p) -> [vtransposelo(in1, in2, 2), vtransposehi(in1, in2, 2)],
    computeType := self >> TVect(TVect(T_Real(32), 2), 2),
));

Class(vextract_half, VecExp_2.unary(), rec(
    computeType := self >> TVect(T_Real(32), 2),
    semantic := (in1, p) -> Checked(Length(in1)=2, in1[p[1]+1]),
    ev := self >> self.t.value(self.semantic(self._vval(1), self.args[2].p)).ev(),
));

Class(vmulcx_half, VecExp_2.binary());
Class(vswapcx_half, VecExp_2.unary());

Class(vload1_half,  VecExp_2.binary());
Class(vstore1_half, VecExpCommand);

# vdup_lane_half(<float32x2>, <int lane_num>) ==> float32x2
Class(vdup_lane_half, VecExp_2.binary(), rec(
    ev := self >> let(lane := self.args[2].ev(),
	Checked(IsInt(lane), lane >=0, lane <= 1,
	        self.t.value( Replicate(2, self.args[1].ev()[1+lane]) )))
));

# vext_half(<float32x2>, <float32x2>, 1) ==> float32x2
# vext_half([a, b], [c, d], 1) = [b c]
Class(vext_half, VecExp_2.ternary(), rec(
    ev := self >> let(num := self.args[3].ev(),
	Checked(IsInt(num), num=1,
	        self.t.value( self.args[1].ev()[2], self.args[2].ev()[1])))
));

#############################################################
# Composite instructions. Serve as a bridge between NEON and
# TL ISA DB. Later rewritten into combination of vuzpq_32f,
# vzipq_32f and vtrnq_32f with vextract_neon_4x32f.
#

Class(vpacklo_neon, _neon_binop_ev, VecExp_4.binary(), rec(
    semantic := (in1, in2, p) -> shuffle(in1, in2, [1, 3, 1, 3], 4, 1)));

Class(vpackhi_neon, _neon_binop_ev, VecExp_4.binary(), rec(
    semantic := (in1, in2, p) -> shuffle(in1, in2, [2, 4, 2, 4], 4, 1)));

Class(vunpacklo_neon, _neon_binop_ev, VecExp_4.binary(), rec(
    semantic := (in1, in2, p) -> unpacklo(in1, in2, 4, 1)));

Class(vunpackhi_neon, _neon_binop_ev, VecExp_4.binary(), rec(
    semantic := (in1, in2, p) -> unpackhi(in1, in2, 4, 1)));

Class(vtransposelo_neon, _neon_binop_ev, VecExp_4.binary(), rec(
    semantic := (in1, in2, p) -> vtransposelo(in1, in2, 4)));

Class(vtransposehi_neon, _neon_binop_ev, VecExp_4.binary(), rec(
    semantic := (in1, in2, p) -> vtransposehi(in1, in2, 4)));

Class(vunpacklolo2_neon, _neon_binop_ev, VecExp_4.binary(), rec(
    semantic := (in1, in2, p) -> unpacklo(in1, in2, 4, 2)));

Class(vunpacklohi2_neon, _neon_binop_ev, VecExp_4.binary(), rec(
    semantic := (in1, in2, p) -> shuffle(in1, in2, [1,2,3,4], 4, 1)));

Class(vunpackhilo2_neon, _neon_binop_ev, VecExp_4.binary(), rec(
    semantic := (in1, in2, p) -> shuffle(in1, in2, [3,4,1,2], 4, 1)));

Class(vunpackhihi2_neon, _neon_binop_ev, VecExp_4.binary(), rec(
    semantic := (in1, in2, p) -> unpackhi(in1, in2, 4, 2)));

Class(vrev_neon, VecExp_4.unary(), rec(
		semantic := (in1, p) -> vrev64(in1,4),
));
#############################################################

Class(vuzpq_32f, _neon_binop_ev, VecExp_4.binary(), rec(
    semantic := (in1, in2, p) -> [shuffle(in1, in2, [1, 3, 1, 3], 4, 1), shuffle(in1, in2, [2, 4, 2, 4], 4, 1)],
    computeType := self >> TVect(TVect(T_Real(32), 4), 2),
));

Class(vzipq_32f, _neon_binop_ev, VecExp_4.binary(), rec(
    semantic := (in1, in2, p) -> [unpacklo(in1, in2, 4, 1), unpackhi(in1, in2, 4, 1)],
    computeType := self >> TVect(TVect(T_Real(32), 4), 2),
));

Class(vtrnq_32f, _neon_binop_ev, VecExp_4.binary(), rec(
    semantic := (in1, in2, p) -> [vtransposelo(in1, in2, 4), vtransposehi(in1, in2, 4)],
    computeType := self >> TVect(TVect(T_Real(32), 4), 2),
));

Class(vextract_neon_4x32f, VecExp_4.unary(), rec(
    computeType := self >> TVect(T_Real(32), 4),
    semantic := (in1, p) -> Checked(Length(in1)=2, in1[p[1]+1]),
    ev := self >> self.t.value(self.semantic(self._vval(1), self.args[2].p)).ev(),
));

Class(vmulcx_neon, VecExp_4.binary());
Class(vswapcx_neon, VecExp_4.unary());

# vload1_neon(<float*>, <float32x4>, <int lane_num>) => float32x4
Class(vload1_neon,  VecExp_4.binary());

# vload_half_neon(<float*>) => float32x2
Class(vload_half_neon, VecExp_2.unary());

# vstore1_neon(<float*>, <float32x4>, <int lane_num>)
Class(vstore1_neon, VecStoreCommand.binary());

# vstore2lo_neon(<float*>, <float32x4>)
Class(vstore2lo_neon, VecStoreCommand.binary());

# vstore2hi_neon(<float*>, <float32x4>)
Class(vstore2hi_neon, VecStoreCommand.binary());

# vcombine_neon(<float32x2>, <float32x>) => float32x4
Class(vcombine_neon, VecExp_4.binary());

# vdup_lane_neon(<float32x4>, <int lane_num>) ==> float32x4
Class(vdup_lane_neon, VecExp_4.binary(), rec(
    ev := self >> let(lane := self.args[2].ev(),
	Checked(IsInt(lane), lane >=0, lane <= 3,
	        self.t.value( Replicate(4, self.args[1].ev()[1+lane]) )))
));

# vext_neon(<float32x4>, <float32x4>, 1) ==> float32x4
# vext_neon([a, b, c, d], [e, f, g, h], 1) = [d e f g]
# vext_neon([a, b, c, d], [e, f, g, h], 2) = [c d e f]
Class(vext_neon, VecExp_4.ternary(), rec(
    ev := self >> let(num := self.args[3].ev(),
	Checked(IsInt(num), num>=0, num<=3,
	        self.t.value( self.args[1].ev(){[5-num .. 4]} :: self.args[2].ev(){[1..4-num]})))
));


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

#F IRDFT(<size>, [<normalized?>]) - Inverse RDFT, non-terminal
#F This is an obsolete inverse real DFT, plese use IPRDFT
#F
Class(IRDFT, NonTerminal, rec(
    abbrevs := [
	N -> Checked(IsInt(N), N > 1, [N, false]),
	(N, norm) -> Checked(IsInt(N), N > 1, IsBool(norm), [N, norm])
    ],
    ndiag := meth(self)
        local Ls, n;
	n := self.params[1];
	if self.params[2] then
	    Ls:=List([1..n], i->2/n);
	    Ls[1]:=1/n;
	    if (n mod 2 = 0) then Ls[n/2+1]:=1/n; fi;
	    return Diag(Ls);
	else return I(self.params[1]);
	fi;
    end,
    transpose := self >> let(exp := 1, Compose(self.ndiag(), RDFT(self.params[1]))),
    terminate := self >> TransposedSPL(TerminateSPL(self.transpose())),
    isReal := True,
    dims := self >> [self.params[1], self.params[1]],
    SmallRandom := () -> Random([2, 3, 4, 6, 8, 12, 16, 24, 32])
));

RulesFor(IRDFT, rec(
    #F IRDFT_Trigonometric: 1984
    #F
    #F   invRDFT_n = blocks * (DCT1_(n/2+1) dirsum DST1_(n/2-1))
    #F
    #F   Wang:
    #F     Fast Algorithms for the Discrete W Transform and the
    #F     Discrete Fourier Transform.
    #F     IEEE Trans. on ASSP, 1984, pp. 803--814
    #F   Britanak/Rao:
    #F     The fast generalized discrete Fourier transforms: A unified
    #F     approach to the discrete sinusoidal transforms computation.
    #F     Signal Processing 79, 1999, pp. 135--150
    #F
    IRDFT_Trigonometric := rec(
	info             := "IRDFT_n -> DCT1_(n/2+1), DST1_(n/2-1)",
	forTransposition := false,
	isApplicable     := P -> P[1] mod 2 = 0 and P[1] > 2,
	allChildren      := P -> [[ DCT1(P[1]/2 + 1), DST1(P[1]/2 - 1) ]],
	rule := (P, C) -> let(
	    n := P[1],
	    exp := 1, # should be an extra IRDFT parameter
	    normalized := P[2],

	    DirectSum(I(1), blocks1(n - 1)) *
	    DirectSum(C[1], exp * C[2] ^ J(n/2 - 1)) *
	    When(normalized,
		Diag(List([0..n-1], c -> When(c mod (n/2) = 0, 1/n, 2/n))),
		I(n)))
    )
));


text:
SymmetricDifference := function(a, b)
  return Union(Difference(a, b), Difference(b, a));
end;

a := ["John", "Serena", "Bob", "Mary", "Serena"];
b := ["Jim", "Mary", "John", "Jim", "Bob"];
SymmetricDifference(a,b);
[ "Jim", "Serena" ]


text:
%Terminals
    DollarSign ::= '$'
    Percent ::= '%'
    _
    a b c d e f g h i j k l m n o p q r s t u v w x y z
%End

%Headers
    /.
        //
        // Each upper case letter is mapped into its corresponding
        // lower case counterpart. For example, if an 'A' appears
        // in the input, it is mapped into %sym_type.%prefix%a%suffix% just
        // like 'a'.
        //
       public  static   tokenKind : number[]=  new Array(128)  ;
         static  __b_init : boolean = %action_type.init_block(%action_type.tokenKind);
        static  init_block(tokenKind : number[]) : boolean
        {
            for (let i = 0; i < tokenKind.length; ++i) {
                tokenKind[i] = 0;
            }
            tokenKind['$'.charCodeAt(0)] = %sym_type.%prefix%DollarSign%suffix%;
            tokenKind['%'.charCodeAt(0)] = %sym_type.%prefix%Percent%suffix%;
            tokenKind['_'.charCodeAt(0)] = %sym_type.%prefix%_%suffix%;

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

            tokenKind['A'.charCodeAt(0)] = %sym_type.%prefix%a%suffix%;
            tokenKind['B'.charCodeAt(0)] = %sym_type.%prefix%b%suffix%;
            tokenKind['C'.charCodeAt(0)] = %sym_type.%prefix%c%suffix%;
            tokenKind['D'.charCodeAt(0)] = %sym_type.%prefix%d%suffix%;
            tokenKind['E'.charCodeAt(0)] = %sym_type.%prefix%e%suffix%;
            tokenKind['F'.charCodeAt(0)] = %sym_type.%prefix%f%suffix%;
            tokenKind['G'.charCodeAt(0)] = %sym_type.%prefix%g%suffix%;
            tokenKind['H'.charCodeAt(0)] = %sym_type.%prefix%h%suffix%;
            tokenKind['I'.charCodeAt(0)] = %sym_type.%prefix%i%suffix%;
            tokenKind['J'.charCodeAt(0)] = %sym_type.%prefix%j%suffix%;
            tokenKind['K'.charCodeAt(0)] = %sym_type.%prefix%k%suffix%;
            tokenKind['L'.charCodeAt(0)] = %sym_type.%prefix%l%suffix%;
            tokenKind['M'.charCodeAt(0)] = %sym_type.%prefix%m%suffix%;
            tokenKind['N'.charCodeAt(0)] = %sym_type.%prefix%n%suffix%;
            tokenKind['O'.charCodeAt(0)] = %sym_type.%prefix%o%suffix%;
            tokenKind['P'.charCodeAt(0)] = %sym_type.%prefix%p%suffix%;
            tokenKind['Q'.charCodeAt(0)] = %sym_type.%prefix%q%suffix%;
            tokenKind['R'.charCodeAt(0)] = %sym_type.%prefix%r%suffix%;
            tokenKind['S'.charCodeAt(0)] = %sym_type.%prefix%s%suffix%;
            tokenKind['T'.charCodeAt(0)] = %sym_type.%prefix%t%suffix%;
            tokenKind['U'.charCodeAt(0)] = %sym_type.%prefix%u%suffix%;
            tokenKind['V'.charCodeAt(0)] = %sym_type.%prefix%v%suffix%;
            tokenKind['W'.charCodeAt(0)] = %sym_type.%prefix%w%suffix%;
            tokenKind['X'.charCodeAt(0)] = %sym_type.%prefix%x%suffix%;
            tokenKind['Y'.charCodeAt(0)] = %sym_type.%prefix%y%suffix%;
            tokenKind['Z'.charCodeAt(0)] = %sym_type.%prefix%z%suffix%;
            return true;
        }

       public  static    getKind(c :number ):number
        {
            return (c < 128 ? %action_type.tokenKind[c] : 0);
        }
    ./
%End


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

DataTypes := rec(
    f64re := rec(
    isFixedPoint := false,
    bits := 64,
    customDataType := "double",
    XType := TPtr(TReal),
    YType := TPtr(TReal),
    TRealCtype := "double",
    includes := [ "<include/omega64.h>" ]
    ),

    f64c := rec(
    isFixedPoint := false,
    bits := 128,
    customDataType := "double_cplx",
    XType := TPtr(TComplex),
    YType := TPtr(TComplex),
    TRealCtype := "double",
    TComplexCtype := "_Complex double",
    includes := [ "<include/omega64c.h>" ]
    ),

    f32re := rec(
    isFixedPoint := false,
    bits := 32,
    customDataType := "float",
    XType := TPtr(TReal),
    YType := TPtr(TReal),
    TRealCtype := "float",
    includes := [ "<include/omega32.h>" ],
    valuePostfix := "f"
    ),

    f32c := rec(
    isFixedPoint := false,
    bits := 64,
    customDataType := "float_cplx",
    XType := TPtr(TComplex),
    YType := TPtr(TComplex),
    TRealCtype := "float",
    TComplexCtype := "_Complex float",
    includes := [ "<include/omega32c.h>" ],
    valuePostfix := "f"
    ),

    i32re := rec(
    isFixedPoint := true,
    bits := 32,
    fracbits := 14,
    customDataType := "int_fp14",
    XType := TPtr(TInt),
    YType := TPtr(TInt),
    TIntCtype := "int",
    TRealCtype := "int",
    includes := [ "<include/omega32i.h>" ]
    ),

    i16re := rec(
    isFixedPoint := true,
    bits := 16,
    fracbits := 7,
    customDataType := "short_fp7",
    XType := TPtr(TInt),
    YType := TPtr(TInt),
    TIntCtype := "signed short",
    TRealCtype := "signed short",
    includes := [ "<include/omega16i.h>" ]
    ),

    i8re := rec(
    isFixedPoint := true,
    bits := 8,
    fracbits := 4,
    customDataType := "char_fp4",
    XType := TPtr(TInt),
    YType := TPtr(TInt),
    TIntCtype := "signed char",
    TRealCtype := "signed char",
    includes := [ "<include/omega8i.h>" ]
    )
);

#F InitDataType(<opts>, <data_type>)
#F   Changes the active data type in <opts> and configures the rest of the
#F   system to generate code with that data type.
#F
#F   data_type = "f64re" | "f32re" | "i32re" | "i16re" | "i8re"
#F
#F   Note: currently this only works correctly with CUnparser, and therefore
#F         best config to use is the one produced by InitIter
#F
#F   Example: opts := InitIter(false);
#F            opts := InitDataType(opts, "i16re");
##
InitDataType := function(opts, data_type)
    local dt;
    Constraint(IsString(data_type));
    dt := DataTypes.(data_type);
    opts := CopyFields(opts, dt);
#    if dt.isFixedPoint then
#   opts.compileStrategy := Concatenation(opts.compileStrategy,
#       [ c -> FixedPointCode(c, dt.bits, dt.fracbits) ]);
#    fi;
    return opts;
end;


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

RulesFor(DFT, rec(
    DFT_SR_Reg := rec(
	isApplicable := P -> P[1] mod 2 = 0 and P[1] > 4,

	allChildren := P -> let(k:=P[2],
	    Map2(Filtered(DivisorPairs(P[1]), dd->dd[2]>2 and (dd[2] mod 2 = 0)),
	    (m,n) -> [DFT(m,k), DFT(2*m,k), UDFT(n,k)])),

	rule := (P,C) -> let( # we know that n >= 4
	    N := P[1], m:=Cols(C[1]), n:=Cols(C[3]), k:=P[2], hf:=(n-2)/2,
	    i := Ind(hf), j := Ind(hf),

	    SUM(        _hhi(N, 2*m, 0,      n/2, C[2] * L(2*m,2)),
		ISum(i, _hhi(N,   m, i+1,    n,   C[1] * Diag(Twid(N,m,k,0,0,i+1)))),
		ISum(j, _hhi(N,   m, j+2+hf, n,   C[1] * Diag(Twid(N,m,k,0,0,j+2+hf))))
	    ),
	    Tensor(I(m), C[3]) *
	    L(N,m)
	)
    ),

    DFT_SR_Reg2 := rec(
	isApplicable := P -> P[1] mod 4 = 0 and P[1] > 4,

	allChildren := P -> let(k:=P[2],
	    Map2(Filtered(DivisorPairs(P[1]), dd->dd[1]>2 and dd[2]>2 and (dd[1] mod 2 = 0) and (dd[2] mod 2 = 0)),
	    (m,n) -> [UDFT(m,k).transpose(), UUDFT(2*m,k).transpose(),
		      UDFT(n,k),             UUDFT(2*n, k)])),

	rule := (P,C) -> let( # we know that n >= 4
	    N := P[1], m:=Cols(C[1]), n:=Cols(C[3]), k:=P[2], nhf:=(n-2)/2, mhf:=(m-2)/2,
	    i := Ind(nhf), j := Ind(nhf), ii := Ind(mhf), jj := Ind(mhf),
	    q := Ind(2*n),
	    dd := Diag(Lambda(q, cond(eq(q,3*n/2), (-E(4))^k, 1))),

	    SUM(        _hhi(N, 2*m, 0,       n/2, C[2] * L(2*m,2)),
		ISum(i, _hhi(N,   m, i+1,     n,   C[1] * Diag(UTwid(N,m,k,0,0,i+1)))),
		ISum(j, _hhi(N,   m, j+2+nhf, n,   C[1] * Diag(UTwid(N,m,k,0,0,j+2+nhf))))
	    ) *
	    SUM(Scat(fStack(H(N,n,0,1), H(N,n,N/2,1))) * dd * C[4] * Gath(H(N, 2*n, 0, m/2)),
		ISum(ii, _hho(N,   n, n+ii*n,         1,  ii+1,     m,   C[3])),
		ISum(jj, _hho(N,   n, 2*n+mhf*n+jj*n, 1,  jj+mhf+2, m,   C[3]))
	    )
	)
    ),

    DFT_SR123_Corr := rec(
	isApplicable := P -> P[1] mod 4 = 0,
	allChildren := P -> let(k:=P[2],
	    Map2(Filtered(DivisorPairs(P[1]), dd -> (dd[1] mod 2=0) and (dd[2] mod 2 = 0)),
	    (m,n) -> [DFT(m,k), DFT3(m,k), DFT(n,k), DFT2(n, k)])),

	rule := function(P,C)
	    local N,m,n,k,nhf,mhf, i,j,ii,jj,q,dd,isum, tw1a, tw1b, tw2a, tw2b;

	    N := P[1]; m:=Cols(C[1]); n:=Cols(C[3]); k:=P[2]; nhf:=(n-2)/2; mhf:=(m-2)/2;
	    i := Ind(nhf); j := Ind(nhf); ii := Ind(mhf); jj := Ind(mhf);
	    q := Ind(n);
	    dd := Diag(Lambda(q, cond(eq(q,n/2), (-E(4))^k, 1)));
	    isum := (var, expr) -> When(var.range=0, [], ISum(var, expr));

	    return
	    SUM(         _hhi(N,   m, 0,       n,   C[1]),
		isum(i,  _hhi(N,   m, i+1,     n,   C[1] * Diag(fPrecompute(UTwid(N,m,k,0,0,i+1))))),
		         _hhi(N,   m, nhf+1,   n,   C[2]),
		isum(j,  _hhi(N,   m, j+nhf+2, n,   C[1] * Diag(fPrecompute(UTwid(N,m,k,0,0,j+2+nhf)))))
	    ) *
	    SUM(         _hho(N,   n, 0,              1,  0,        m,    C[3]),
		isum(ii, _hho(N,   n, n+ii*n,         1,  ii+1,     m,    C[3])),
		         _hho(N,   n, n+mhf*n,        1,  mhf+1,    m, dd*C[4]),
		isum(jj, _hho(N,   n, 2*n+mhf*n+jj*n, 1,  jj+mhf+2, m,    C[3]))
	    );
	end
    ),

    DFT_SR123 := rec(
	isApplicable := P -> P[1] mod 4 = 0,
	maxRadix := 4,
	allChildren := (self, P) >> let(k:=P[2],
	    Map2(Filtered(DivisorPairs(P[1]), dd -> (dd[1] mod 2=0 and dd[1]<=self.maxRadix) and (dd[2] mod 2 = 0)),
	    (m,n) -> [DFT(m,k), DFT3(m,k), DFT(n,k), DFT2(n, k)])),

	inplace := Inplace,
	rule := meth(self,P,C)
	    local N,m,n,k,nhf,mhf, i,j,ii,jj,q,dd,isum, tw1a, tw1b, tw2a, tw2b, d1;

	    N := P[1]; m:=Cols(C[1]); n:=Cols(C[3]); k:=P[2]; nhf:=(n-2)/2; mhf:=(m-2)/2;
	    i := Ind(nhf); j := Ind(nhf); ii := Ind(mhf); jj := Ind(mhf);
	    q := Ind(n);
	    dd := Diag(Lambda(q, cond(eq(q,n/2), (-E(4))^k, 1)));
	    isum := (var, expr) -> When(var.range=0, [], ISum(var, expr));

	    tw1a := fPrecompute(fCompose(dOmega(N,k), dLin(mhf, i+1, i+1,     TInt)));
	    tw1b := fPrecompute(fCompose(dOmega(N,k), dLin(mhf, i+1, (i+1)*(mhf+2), TInt)));

	    tw2a := fPrecompute(fCompose(dOmega(N,k), dLin(mhf, j+2+nhf, j+2+nhf,     TInt)));
	    tw2b := fPrecompute(fCompose(dOmega(N,k), dLin(mhf, j+2+nhf, (j+2+nhf)*(mhf+2), TInt)));

	    d1 := fConst(1, 1.0);

	    return
	    self.inplace(
	    SUM(         _hhi(N,   m, 0,       n,   C[1]),
		isum(i,  _hhi(N,   m, i+1,     n,   C[1] * When(tw1a.domain()=0, I(2), Diag(diagDirsum(d1, tw1a, d1, tw1b))))),
		         _hhi(N,   m, nhf+1,   n,   C[2]),
		isum(j,  _hhi(N,   m, j+nhf+2, n,   C[1] * When(tw2a.domain()=0, I(2), Diag(diagDirsum(d1, tw2a, d1, tw2b)))))
	    )) *
	    SUM(         _hho(N,   n, 0,              1,  0,        m,    C[3]),
		isum(ii, _hho(N,   n, n+ii*n,         1,  ii+1,     m,    C[3])),
		         _hho(N,   n, n+mhf*n,        1,  mhf+1,    m,    C[4]), #dd * C[4]),
		isum(jj, _hho(N,   n, 2*n+mhf*n+jj*n, 1,  jj+mhf+2, m,    C[3]))
	    );
	end
    )
));

 RulesFor(DFT3, rec(
    DFT3_CT_Twids := Inherit(DFT_CT, rec(
 	switch := false,
	maxRadix := 4,
 	allChildren  := (self, P) >> Map2(Filtered(DivisorPairs(P[1]), dd->dd[1] <= self.maxRadix),
 	    (m,n) -> [ DFT(m, P[2]), DFT3(n, P[2]) ]),

	inplace := Inplace,
 	rule := meth(self, P,C)
	    local N, m, n, k, i, tw;
	    N := P[1];
	    m := Cols(C[1]); n := Cols(C[2]); k := P[2];
	    i := Ind(n);
	    tw := fPrecompute(fCompose(dOmega(2*N,k), dLin(m-1, 2*i+1, 2*i+1,  TInt)));

	    return self.inplace(ISum(i, _hhi(N,  m,  i,  n,  C[1] * Diag(diagDirsum(fConst(1,1.0), tw))))) *
 	           Tensor(I(m), C[2]) *
		   L(N, m);
 	end
     )),

    DFT3_SR34_Corr := rec(
	isApplicable := P -> P[1] mod 4 = 0,
	allChildren := P -> let(k:=P[2],
	    Map2(Filtered(DivisorPairs(P[1]), dd -> (dd[1] mod 2=0) and (dd[2] mod 2 = 0)),
	    (m,n) -> [DFT(m,k), DFT3(n,k), DFT4(n, k)])),

	rule := function(P,C)
	    local N,m,n,k,nhf,mhf, i,j,ii,jj,q,dd,isum, tw1a, tw1b, tw2a, tw2b;

	    N := P[1]; m:=Cols(C[1]); n:=Cols(C[3]); k:=P[2]; nhf:=(n-2)/2; mhf:=(m-2)/2;
	    i := Ind(n); ii := Ind(mhf); jj := Ind(mhf);
	    isum := (var, expr) -> When(var.range=0, [], ISum(var, expr));

	    return
            isum(i,  _hhi(N,   m, i,     n,   C[1] * Diag(fPrecompute(UTwid1(N,m,k,1/2,0,i))))) *
	    SUM(         _hho(N,   n, 0,              1,  0,        m,    C[2]),
		isum(ii, _hho(N,   n, n+ii*n,         1,  ii+1,     m,    C[2])),
		         _hho(N,   n, n+mhf*n,        1,  mhf+1,    m,    C[3]),
		isum(jj, _hho(N,   n, 2*n+mhf*n+jj*n, 1,  jj+mhf+2, m,    C[2]))
	    );
	end
    ),

 ));

RulesFor(UUDFT, rec(
   UUDFT1_Base4 := BaseRule(UUDFT1, [4, @])
));

RulesFor(UDFT1, rec(
   UDFT1_Base2 := rec(
       isApplicable := P -> P[1] = 2,
       rule := (P, C) -> I(2)
   ),

   UDFT1_Base4 := rec(
       isApplicable := P -> P[1] = 4,
       rule := (P, C) ->
           DirectSum(I(2), F(2))^L(4,2) *
	   Diag(Tw1(4, 2, 1)) *
	   Tensor(I(2), F(2)) *
	   L(4, 2)
   ),

   UDFT1_toDFT := rec(
       isApplicable := P -> P[1] > 4,
       allChildren := P -> [[ DFT(P[1], P[2]) ]],
       rule := (P, C) ->
           DirectSum((1/2)*F(2), I(P[1]-2))^L(P[1],P[1]/2) * C[1]
   ),

   UDFT1_SR_Reg := rec(
	isApplicable := P -> Is2Power(P[1]) and P[1] > 4,

	allChildren := P -> let(k:=P[2], Map2(Filtered(DivisorPairs(P[1]), dd->dd[2]>2),
	    (m,n) -> [DFT(m,k), UDFT(2*m,k), UDFT(n,k)])),

	forTransposition := false,

	rule := (P,C) -> let( # we know that n >= 4
	    N := P[1], m:=Cols(C[1]), n:=Cols(C[3]), k:=P[2], j := Ind((n-2)/2), jj := Ind((n-2)/2),
	    SUM(
		Scat(H(N,2*m,0,n/2)) * C[2] * L(2*m,2) * Gath(H(N,2*m,0,n/2)),
		ISum(j,
		     Scat(H(N,m,j+1,n)) *
		     C[1] * Diag(Twid(N,m,k,0,0,j+1)) *
		     Gath(H(N, m,j+1,n))),
		ISum(jj,
		     Scat(H(N,m,jj+2+(n-2)/2,n)) *
		     C[1] * Diag(Twid(N,m,k,0,0,jj+2+(n-2)/2)) *
		     Gath(H(N, m,jj+2+(n-2)/2,n)))

	    ) *
	    Tensor(I(m), C[3]) * L(N,m))
    )
));

#Unbind(UDFT.transpose);
#SwitchRules(DFT, [1,3,18]);
#DFT_CT.isApplicable := P -> P[1]=4;


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

# This is the code generator I use for generating HDL descriptions.  The
# output from this is sent to my back-end.

Class(HDLCodegen, Codegen, rec(
    initXY := function(x, y, opts)
        if IsBound(opts.XType) and not IsList(x) then
            if not IsArrayT(opts.XType) and not IsPtrT(opts.XType)
                then Error("opts.XType must be a pointer or array type. This has recently changed.",
                           "If you used TReal before, use TPtr(TReal) now"); fi;
            x.t := opts.XType;
            x.t := When(IsBound(opts.useRestrict) and opts.useRestrict, x.t.restrict(), x.t);
        fi;
        if IsBound(opts.YType) and not IsList(y) then
            if not IsArrayT(opts.YType) and not IsPtrT(opts.YType)
                then Error("opts.YType must be a pointer or array type. This has recently changed.",
                           "If you used TReal before, use TPtr(TReal) now"); fi;
            y.t := opts.YType;
            y.t := When(IsBound(opts.useRestrict) and opts.useRestrict, y.t.restrict(), y.t);
        fi;
        return [x, y];
    end,

    Formula := meth(self, o, y, x, opts)

        local code, datas, prog, params, sub, initsub;

	[x, y] := self.initXY(x,y,opts);

	o := o.child(1);
	params := Set(Collect(o, param));

	datas := Collect(o, FDataOfs);
	o := BlockSums(opts.globalUnrolling, o);
	code := self(o, y, x, opts);
    code := ESReduce(code, opts);
	code := RemoveAssignAcc(code);
	code := BlockUnroll(code, opts);
        # code := PowerOpt(code);
	code := DeclareHidden(code);
	if IsBound(opts.isFixedPoint) and opts.isFixedPoint then
	    code := FixedPointCode(code, opts.bits, opts.fracbits);
	fi;

	sub := Cond(IsBound(opts.subName), opts.subName, "transform");
	initsub := Cond(IsBound(opts.subName), Concat("init_", opts.subName), "init");
#	code := func(TVoid, sub, Concatenation(params, [y, x]), code);
#	code :=

#	if IsBound(opts.generateInitFunc) and opts.generateInitFunc then
#	    prog := program(
#		decl(List(datas, x->x.var),
#		    chain(
#			func(TVoid, initsub, [], chain(List(datas, x -> SReduce(x.var.init, opts)))),
#			code
#		    )));
#	else
#	    prog := program(code);
#	fi;
	prog := code;
	prog.dimensions := o.dimensions;
	return prog;
    end,

#    BB := (self,o,y,x,opts) >> MarkForUnrolling(self(o.child(1), y, x, opts)),
    BB := (self,o,y,x,opts) >> MarkForUnrolling(
        When(IsBound(o.bbnum),
            o.bbnum,
            0
        ),
        self(o.child(1), y, x, opts)
    ),

    Buf := (self,o,y,x,opts) >> self(o.child(1), y, x, opts),
    Grp := (self,o,y,x,opts) >> self(o.child(1), y, x, opts),

    COND := (self,o,y,x,opts) >> IF(
        When(IsFunction(o.cond), o.cond.at(0), o.cond),
        self(o.child(1), y, x, opts), self(o.child(2), y, x, opts)),

#    COND := (self,o,y,x,opts) >> IF(o.cond, self(o.child(1), y, x, opts), self(o.child(2), y, x, opts)),

    Diag := (self, o, y, x, opts) >> let(i := Ind(), elt := o.element.lambda(),
	loop(i, elt.domain(), assign(nth(y,i), elt.at(i) * nth(x,i)))),

    RCDiag := (self, o, y, x, opts) >> let(i := Ind(), elt := o.element.lambda(),
	re := elt.at(2*i), im := elt.at(2*i+1),
	loop(i, elt.domain()/2, chain(
		assign(nth(y,2*i),   re * nth(x,2*i) - im * nth(x,2*i+1)),
		assign(nth(y,2*i+1), im * nth(x,2*i) + re * nth(x,2*i+1))))),

    ColVec := (self, o, y, x, opts) >> let(i := Ind(), func := o.element.lambda(),
	loop(i, func.domain(), assign(nth(y,i), mul(func.at(i), nth(x,0))))),

    RowVec := (self, o, y, x, opts) >> let(i := Ind(), func := o.element.lambda(),
	t := TempVar(x.t.t),
	chain(assign(t,0),
	    loop(i, func.domain(), assign(t, add(t, mul(func.at(i), nth(x,i))))),
	    assign(nth(y,0), t))),

    Scale := (self, o, y, x, opts) >> let(i := Ind(),
	chain(self(o.child(1), y, x, opts),
	    loop(i, Rows(o), assign(nth(y,i), mul(o.scalar, nth(y,i)))))),

    I := (self, o, y, x, opts) >> let(i := Ind(Rows(o)),
	loop(i, i.range, assign(nth(y,i), nth(x,i)))),

    Blk := (self, o, y, x, opts) >> Cond(
	Rows(o)=2 and Cols(o)=2, Blk2code(o, y, x),
        Rows(o)=4 and Cols(o)=4, Blk4code(o, y, x),
	let(j:=Ind(), i:=Ind(), t:=TempVar(x.t.t), mat:=V(o.element), d:=Dat(mat.t),
	    data(d, mat,
		loop(j, Rows(o), decl(t, chain(
		     assign(t, 0),
		     loop(i, Cols(o), assign(t, t + nth(nth(d,j),i) * nth(x,i))),
		     assign(nth(y,j), t))))))),

    toeplitz := (self, o, y, x, opts) >> self.Blk(o.obj, y, x, opts),

    Blk1 := (self, o, y, x, opts) >> assign(nth(y,0), mul(toExpArg(o.element), nth(x,0))),

    BlkConj := (self, o, y, x, opts) >> assign(nth(y,0), conj(nth(x,0))),

    Prm := (self, o, y, x, opts) >> let(i:=Ind(), func:=o.func.lambda(),
	loop(i, Rows(o), assign(nth(y, i), nth(x, func.at(i))))),

    L := meth(self, o, y, x, opts)
        return let(i:=Ind(), func:=o.lambda(),
            loop(i, Rows(o), assign(nth(y, i), nth(x, func.at(i)))).unroll());
    end,

    O := (self, o, y, x, opts) >> let(i:=Ind(),
	loop(i, o.params[1], assign(nth(y, i), V(0)))),

    DirectSum := (self,o,y,x,opts) >> self(o.sums(), y, x, opts),
    Tensor := (self,o,y,x,opts) >> self(o.sums(), y, x, opts),
    IDirSum := (self,o,y,x,opts) >> self(o.sums(), y, x, opts),

    Gath := meth(self, o, y, x, opts)
        local i, func, rfunc, ix;
	i := Ind(); func := o.func.lambda();

	if IsBound(o.func.rlambda) then
	    rfunc := o.func.rlambda();
	    ix := var.fresh_t("ix", TInt);
	    return decl(ix, chain(
		assign(ix, func.at(0)),
		assign(nth(y, 0), nth(x, ix)),
		loop(i, o.func.domain()-1,
		    chain(assign(ix, rfunc.at(ix)),
			assign(nth(y, i+1), nth(x, ix))))));
	else
	    return loop(i, o.func.domain(), assign(nth(y,i), nth(x, func.at(i))));
	fi;
    end,

    Scat := meth(self, o, y, x, opts)
        local i, func, rfunc, ix;
	i := Ind(); func := o.func.lambda();
	if IsBound(o.func.rlambda) then
	    rfunc := o.func.rlambda();
	    ix := var.fresh_t("ix", TInt);
	    return decl(ix, chain(
		assign(ix, func.at(0)),
		assign(nth(y, ix), nth(x, 0)),
		loop(i, o.func.domain()-1,
		    chain(assign(ix, rfunc.at(ix)),
			assign(nth(y, ix), nth(x, i+1))))));
	else
	    return loop(i, o.func.domain(), assign(nth(y,func.at(i)), nth(x, i)));
	fi;
    end,

    SUM := (self, o, y, x, opts) >> chain(List(o.children(), c -> self(c, y, x, opts))),

    ISum := (self, o, y, x, opts) >> let(myloop := When(IsSymbolic(o.domain), loopn, loop),
    myloop(o.var, o.domain,
        self(o.child(1), y, x, opts))),

    # NOTE: get rid of _acc
    SUMAcc := (self, o, y, x, opts) >> let(ii := Ind(),
       When(not Same(ObjId(o.child(1)), Gath),  # NOTE: come up with a general condition
       chain(
           loop(ii, Rows(o), assign(nth(y, ii), V(0))),
           List(o.children(), c -> _acc(self(c, y, x, opts), y))),
       chain(
           self(o.child(1), y, x, opts),
           List(Drop(o.children(), 1), c -> _acc(self(c, y, x, opts), y))))),

    ISumAcc := (self, o, y, x, opts) >> let(ii := Ind(), chain(
    loop(ii, Rows(o), assign(nth(y, ii), V(0))),
    loop(o.var, o.domain, _acc(self(o.child(1), y, x, opts), y)))),

#     Compose := meth(self, o,y,x,opts)
#         local ch, numch, vecs, allow, i;
# #   if IsBound(opts._inplace) and opts._inplace then
# #       return chain(
# #       List(Reversed(o.children()), c -> self(c, x, x, opts))); fi;
#         #VIENNA filtering DMPGath, DMPScat out here because they do not generate code
#     ch := Filtered(o.children(), i-> not i.name in ["DMPGath","DMPScat"]);
#     numch := Length(ch);
#     vecs := [y];
#     allow := (x<>y);
#     for i in [1..numch-1] do
#             if allow and ObjId(ch[i])=Inplace then vecs[i+1] := vecs[i];
#         else vecs[i+1] := TempVec(TArray(TempArrayType(y, x), Cols(ch[i])));
#         fi;
#     od;
#     vecs[numch+1] := x;
#     for i in Reversed([1..numch]) do
#             if allow and ObjId(ch[i])=Inplace
#         then vecs[i] := vecs[i+1]; fi;
#     od;

#         # everything was inplace, make it go from x -> y as expected
#     if vecs[1] = vecs[numch+1] then vecs[1] := y; fi;
#     if vecs[1] = vecs[numch+1] then vecs[numch+1] := x; fi;

#     [vecs, ch] := [Reversed(vecs), Reversed(ch)];
#     return decl( Difference(vecs{[2..Length(vecs)-1]}, [x,y]),
#         chain( List([1..numch], i -> When(vecs[i+1]=vecs[i],
#             self(ch[i], vecs[i],   vecs[i], CopyFields(opts, rec(_inplace:=true))),
#             self(ch[i], vecs[i+1], vecs[i], opts)))));
#     end,

    Compose := meth(self, o,y,x,opts)
        local ch, numch, vecs, i, cmd, code;
        ch    := o.children();
        numch := Length(ch);
        vecs  := [y];

        # we like to evaluate right (last) to left (first), where the values
        # in parens refer to the position of the elements in the array.
        #
        # So here, we start at the output (first entry in array) and
        # walk towards the input, creating temporary arrays as necessary
        for i in [1..numch-1] do
            vecs[i+1] := When(ObjId(ch[i]) = Inplace,
                vecs[i],
                TempArray(y,x,ch[i])
            );
        od;

        # the last entry must be the input.
        vecs[numch+1] := x;

        # now we walk in the opposite direction, copying through
        # input as far as we can.
        for i in Reversed([1..numch]) do
            if ObjId(ch[i]) = Inplace then
                vecs[i] := vecs[i+1];
            fi;
        od;

        # If all children were evaluated inplace, the output will be in x
        # Make it go from x -> y as expected
        if vecs[1] = vecs[numch+1] then vecs[1] := y; fi;
        if vecs[1] = vecs[numch+1] then vecs[numch+1] := x; fi;

        # order them so that first to be evaluated is first in array.
        vecs := Reversed(vecs);
        ch   := Reversed(ch);

        # Wrap code in variable declaration. Each entry in vecs will contain multiple
        # arrays in the case of multi-input/output (i.e. OL)
        return decl(
            Difference(Flat(vecs{[2..Length(vecs)-1]}), Flat([x,y])),
            chain(
                List([1..numch], i ->
                     self(ch[i], vecs[i+1], vecs[i], opts)))
        );
    end,

    Inplace := (self, o, y, x, opts) >>
#        self(o.child(1), x, x, CopyFields(opts, rec(_inplace:=true))),
        self(o.child(1), y, x, opts), # Compose will handle these somehow

    Data := meth(self, o, y, x, opts)
        local val;
        o.var.isData := true;
        val := When(IsFunction(o.value), o.value.tolist(), o.value);
        val := When(IsValue(val), val, o.var.t.value(val));
        return data(o.var, val, self(o.child(1), y, x, opts));
    end,

#     Data := (self, o, y, x, opts) >> decl(o.var, chain(
# 	    assign(o.var, o.value.tolist()),
# 	    self(o.child(1), y, x, opts))
#     ),

    #   NOTE: use pointers to do pingpong?
    ICompose := (self, o, y, x, opts) >>
        chain([let(
            t :=  Dat1d(x.t.t, Rows(o)),
            newind := Ind(o.domain),
            decl([t], chain(
            When(IsOddInt(o.domain),
            SubstVars(Copy(self(o.child(1), y, x, opts)), tab((o.var.id) := (V(0)))),
                skip()
            ),
            When(IsEvenInt(o.domain),
                chain(
                    SubstVars(Copy(self(o.child(1), t, x, opts)), tab((o.var.id) := (V(0)))),
                    SubstVars(Copy(self(o.child(1), y, t, opts)), tab((o.var.id) := (V(1))))
                ),
                skip()
            ),
            loop(newind, (o.domain/2),
                chain([
                    SubstVars(Copy(self(o.child(1), x, y, opts)), tab((o.var.id) := (2*newind)+2)),
                    SubstVars(Copy(self(o.child(1), y, x, opts)), tab((o.var.id) := (2*newind)+3))
                ])
            )
            )))])
));


text:
IsDeranged := function(a, b)
	local i, n;
	for i in [1 .. Size(a)] do
		if a[i] = b[i] then
			return false;
		fi;
	od;
	return true;
end;

# This solution will find all deranged pairs of any length.
Deranged := function(name)
	local sol, ana, u, v;
	sol := [ ];
	ana := Anagrams(name);
	for u in ana do
		for v in Combinations(u, 2) do
			if IsDeranged(v[1], v[2]) then
				Add(sol, v);
			fi;
		od;
	od;
	return sol;
end;

# Now we find all deranged pairs of maximum length
a := Deranged("unixdict.txt");;
n := Maximum(List(a, x -> Size(x[1])));
Filtered(a, x -> Size(x[1]) = n);
# [ [ "excitation", "intoxicate" ] ]


text:
# This is another ugly hack to make the GAP Help System
# play ball. Let us please fix this soon.
# TODO: This is now broken because we got rid of parsing
#       on the python side. HELP now should result
#       in a record that can be sent back to jupyter
#       as a JSON string
HELP_VIEWER_INFO.jupyter_online :=
    rec(
         type := "url",
         show := function( data )
             # data[1] is the text preceding the hyperlink (name of the help book),
             # data[2] is the text to be linked, and data[3] is the URL
             local p,r;

             p := data[3];

             for r in GAPInfo.RootPaths do
                 p := ReplacedString(data[3], r, "https://www.gap-system.org/Manuals/");
             od;
             return JupyterRenderable( rec( ("text/html") := Concatenation( data[1], ": ", data[2], "") )
                                     , rec( ) );
         end
        );

HELP_VIEWER_INFO.jupyter_local :=
    rec( type := "url",
         show := function( data )
             # data[1] is the text preceding the hyperlink (name of the help book),
             # data[2] is the text to be linked, and data[3] is the URL
             local p,r;

             p := data[3];

             for r in GAPInfo.RootPaths do
                 p := ReplacedString(data[3], r, "/");
             od;
             return JupyterRenderable( rec( ("text/html") := Concatenation( data[1], ": ", data[2], "") )
                                     , rec( ) );
         end);

#############################################################################
##
#F  GET_HELP_URL( <match> ) . . . . . .  print the url for the help section
##
##  Based on HELP_PRINT_MATCH
##
##  <match> is [book, entrynr]
##
InstallGlobalFunction(GET_HELP_URL, function(match)
local book, entrynr, viewer, hv, pos, type, data;
  book := HELP_BOOK_INFO(match[1]);
  entrynr := match[2];
  viewer:= UserPreference("HelpViewers");
  if HELP_LAST.NEXT_VIEWER = false then
    hv := viewer;
  else
    pos := Position( viewer, HELP_LAST.VIEWER );
    if pos = fail then
      hv := viewer;
    else
      hv := viewer{Concatenation([pos+1..Length(viewer)],[1..pos])};
    fi;
    HELP_LAST.NEXT_VIEWER := false;
  fi;
  for viewer in hv do
    # type of data we need now depends on help viewer
    type := HELP_VIEWER_INFO.(viewer).type;
    # get the data via appropriate handler
    data := HELP_BOOK_HANDLER.(book.handler).HelpData(book, entrynr, type);
    if data <> fail then
      # show the data
      return HELP_VIEWER_INFO.(viewer).show(
        [ book.bookname, StripEscapeSequences(book.entries[entrynr][1]), data]);
          # name of the help book, the text to be linked, and the URL
    else
        return JupyterRenderable( rec( ("text/html") := Concatenation( book.bookname, ": "
                                                                       , StripEscapeSequences(book.entries[entrynr][1])
                                                                       , " - no html help available. Please check other formats!" ) )
                                , rec( ) );
    fi;
    HELP_LAST.VIEWER := viewer;
  od;
  HELP_LAST.BOOK := book;
  HELP_LAST.MATCH := entrynr;
  HELP_LAST.VIEWER := viewer;
  return true;
end);

InstallGlobalFunction(JUPYTER_HELP_SHOW_MATCHES, function( books, topic, frombegin )
local   exact,  match,  x,  lines,  cnt,  i,  str,  n, res;

  # first get lists of exact and other matches
  x := HELP_GET_MATCHES( books, topic, frombegin );
  exact := x[1];
  match := x[2];

  # no topic found
  if 0 = Length(match) and 0 = Length(exact)  then
    Print( "Help: no matching entry found\n" );
    return false;

  # one exact or together one topic found
  elif 1 = Length(exact) or (0 = Length(exact) and 1 = Length(match)) then
    if Length(exact) = 0 then exact := match; fi;
    i := exact[1];
    return GET_HELP_URL(i);

  # more than one topic found, show overview in pager
  else
    lines :=
        ["","Help: several entries match this topic - type ?2 to get match [2]\n"];
        # there is an empty line in the beginning since `tail' will start from line 2
    HELP_LAST.TOPICS:=[];
    cnt := 0;
    # show exact matches first
    match := Concatenation(exact, match);
    res:="";
    for i  in match  do
      cnt := cnt+1;
      topic := Concatenation(i[1].bookname,": ",i[1].entries[i[2]][1]);
      Add(HELP_LAST.TOPICS, i);
      Append(res, GET_HELP_URL(i)!.data.("text/html"));
      Append(res, "
");
    od;
    return JupyterRenderable( rec( ("text/html") := res )
                            , rec( ) );
  fi;
end);

InstallGlobalFunction(JUPYTER_HELP, function( str )
  local origstr, nwostr, p, book, books, move, add;

  origstr := ShallowCopy(str);
  nwostr := NormalizedWhitespace(origstr);

  # extract the book
  p := Position( str, ':' );
  if p <> fail  then
      book := str{[1..p-1]};
      str  := str{[p+1..Length(str)]};
  else
      book := "";
  fi;

  # normalizing for search
  book := SIMPLE_STRING(book);
  str := SIMPLE_STRING(str);

  # we check if `book' MATCH_BEGINs some of the available books
  books := Filtered(HELP_KNOWN_BOOKS[1], bn-> MATCH_BEGIN(bn, book));
  if Length(book) > 0 and Length(books) = 0 then
    Print("Help: None of the available books matches (try: '?books').\n");
    return;
  fi;

  # function to add a topic to the ring
  move := false;
  add  := function( books, topic )
      if not move  then
          HELP_RING_IDX := (HELP_RING_IDX+1) mod HELP_RING_SIZE;
          HELP_BOOK_RING[HELP_RING_IDX+1]  := books;
          HELP_TOPIC_RING[HELP_RING_IDX+1] := topic;
      fi;
  end;

  # if the topic is empty show the last shown one again
  if  book = "" and str = ""  then
       if HELP_LAST.BOOK = 0 then
         HELP("Tutorial: Help");
       else
         return GET_HELP_URL( [HELP_LAST.BOOK, HELP_LAST.MATCH] );
       fi;
       return;

  # if topic is "&" shobn;w last topic again, but with next viewer in viewer
  # list, or with last viewer again if there is no next one
  elif book = "" and str = "&" and Length(nwostr) = 1 then
       if HELP_LAST.BOOK = 0 then
         HELP("Tutorial: Help");
       else
         HELP_LAST.NEXT_VIEWER := true;
         return GET_HELP_URL( [HELP_LAST.BOOK, HELP_LAST.MATCH] );
       fi;
       return;

  # if the topic is '-' we are interested in the previous search again
  elif book = "" and str = "-" and Length(nwostr) = 1  then
      HELP_RING_IDX := (HELP_RING_IDX-1) mod HELP_RING_SIZE;
      books := HELP_BOOK_RING[HELP_RING_IDX+1];
      str  := HELP_TOPIC_RING[HELP_RING_IDX+1];
      move := true;

  # if the topic is '+' we are interested in the last section again
  elif book = "" and str = "+" and Length(nwostr) = 1  then
      HELP_RING_IDX := (HELP_RING_IDX+1) mod HELP_RING_SIZE;
      books := HELP_BOOK_RING[HELP_RING_IDX+1];
      str  := HELP_TOPIC_RING[HELP_RING_IDX+1];
      move := true;
  fi;

  # number means topic from HELP_LAST.TOPICS list
  if book = "" and ForAll(str, a-> a in "0123456789") then
      HELP_SHOW_FROM_LAST_TOPICS(Int(str));

  # if the topic is '<' we are interested in the one before 'LastTopic'
  elif book = "" and str = "<" and Length(nwostr) = 1  then
      HELP_SHOW_PREV();

  # if the topic is '>' we are interested in the one after 'LastTopic'
  elif book = "" and str = ">" and Length(nwostr) = 1  then
      HELP_SHOW_NEXT();

  # if the topic is '<<' we are interested in the previous chapter intro
  elif book = "" and str = "<<"  then
      HELP_SHOW_PREV_CHAPTER();

  # if the topic is '>>' we are interested in the next chapter intro
  elif book = "" and str = ">>"  then
      HELP_SHOW_NEXT_CHAPTER();

  # if the subject is 'Welcome to GAP' display a welcome message
  elif book = "" and str = "welcome to gap"  then
      if HELP_SHOW_WELCOME(book)  then
          add( books, "Welcome to GAP" );
      fi;

  # if the topic is 'books' display the table of books
  elif book = "" and str = "books"  then
      if HELP_SHOW_BOOKS()  then
          add( books, "books" );
      fi;

  # if the topic is 'chapters' display the table of chapters
  elif str = "chapters"  or str = "contents" or book <> "" and str = "" then
      if ForAll(books, b->  HELP_SHOW_CHAPTERS(b)) then
        add( books, "chapters" );
      fi;

  # if the topic is 'sections' display the table of sections
  elif str = "sections"  then
      if ForAll(books, b-> HELP_SHOW_SECTIONS(b)) then
        add(books, "sections");
      fi;

  # if the topic is '?<string>' search the index for any entries for
  # which <string> is a substring (as opposed to an abbreviation)
  elif Length(str) > 0 and str[1] = '?'  then
      str := str{[2..Length(str)]};
      NormalizeWhitespace(str);
      return HELP_SHOW_MATCHES( books, str, false);

  # search for this topic
  elif IsJupyterRenderable( HELP_SHOW_MATCHES( books, str, true ) ) then
      return HELP_SHOW_MATCHES( books, str, true );
  elif origstr in NAMES_SYSTEM_GVARS then
      Print( "Help: '", origstr, "' is currently undocumented.\n",
             "      For details, try ?Undocumented Variables\n" );
  elif book = "" and
                 ForAny(HELP_KNOWN_BOOKS[1], bk -> MATCH_BEGIN(bk, str)) then
      Print( "Help: Are you looking for a certain book? (Trying '?", origstr,
             ":' ...\n");
      return HELP( Concatenation(origstr, ":") );
  else
     # seems unnecessary, since some message is already printed in all
     # cases above (?):
     # Print( "Help: Sorry, could not find a match for '", origstr, "'.\n");
  fi;
end);

# Load some help stuff (Experimental)
InstallGlobalFunction(JUPYTER_FindManSection,
function(file, name)
    local xml, sections, p, s, res;
    xml := ParseTreeXMLFile(file);
    CheckAndCleanGapDocTree(xml);
    sections := XMLElements(xml, "ManSection");;
    res := [];
    for s in sections do
        if IsBound(s.content) then
            p := PositionProperty(s.content, x ->
                                               IsBound(x.attributes) and
                                             IsBound(x.attributes.Name) and
                                             x.attributes.Name = name);
            if p <> fail then
                Add(res, s);
            fi;
        fi;
    od;
    return res;
end);

max_stars_repo_name: ZachNewbery/JupyterKernel

max_issues_repo_name: ZachNewbery/JupyterKernel

max_forks_repo_name: ZachNewbery/JupyterKernel


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

#######################################################################################
#   AltiVec 4-way 32-bit float instructions

# shuffle operations
Class(vbinop_4x32f_av, vbinop_av, rec(
    v := self >> 4,
    #ctype := self >> "float",
)); #Class added for AltiVec

Class(vunpacklo_4x32f_av, vbinop_4x32f_av, rec(
    altivec := "vec_mergeh",
    semantic := (in1, in2, p) -> unpacklo(in1, in2, 4, 1)
)); #Class added for AltiVec

Class(vunpackhi_4x32f_av, vbinop_4x32f_av, rec(
    altivec := "vec_mergel",
    semantic := (in1, in2, p) -> unpackhi(in1, in2, 4, 1)
)); #Class added for AltiVec

#   binary instructions
vunpacklo_4x32f.altivec := "vec_mergel";
vunpackhi_4x32f.altivec := "vec_mergeh";

Class(vperm_4x32f, vbinop_4x32f_av, rec(
    altivec := "vec_perm",
    semantic := (in1, in2, p) -> vpermop(in1, in2, p, 4),
    params := self >> sparams(4,8),
    permparams := aperm
));

#   unary instructions
Class(vuperm_4x32f, vunbinop_av, rec(
    binop := vperm_4x32f
));

#######################################################################################
#   AltiVec 8-way 16-bit integer instructions

#   binary instructions
vunpacklo_8x16i.altivec := "vec_mergel";
vunpackhi_8x16i.altivec := "vec_mergeh";

Class(vperm_8x16i, VecExp_8.binary(), rec(
    altivec := "vec_perm",
    semantic := (in1, in2, p) -> vpermop(in1, in2, p, 8),
    params := self >> sparams(8,16),
    permparams := aperm
));

#   unary instructions
Class(vuperm_8x16i, VecExp_8.unaryFromBinop(vperm_8x16i));

#######################################################################################
#   AltiVec 16-way 8-bit integer instructions

#   binary instructions
vunpacklo_16x8i.altivec := "vec_mergel";
vunpackhi_16x8i.altivec := "vec_mergeh";

Class(vperm_16x8i, VecExp_8.binary(), rec(
    altivec := "vec_perm",
    semantic := (in1, in2, p) -> vpermop(in1, in2, p, 16),
    params := self >> sparams(16,32),
    permparams := aperm
));

#   unary instructions
Class(vuperm_16x8i, VecExp_8.unaryFromBinop(vperm_16x8i));


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

TempArrayType := (child, y, x, index) ->
let(When(IsBound(child.a.t_in), child.a.t_in[index],
		let( X := Flat([x])[1], Y:= Flat([y])[1],
			Cond(
				IsBound(X.t.t) and X.t.t = TComplex, 		TComplex,
				IsBound(Y.t.t) and Y.t.t = TComplex, 		TComplex,
				ObjId(child) = ISumReduction, 				child.idval.t,
				ObjId(child) = PointWise, 					child.op.expr.t,
				X.t.t))));
TempArray := (y, x, child) -> let(cols := Flat([ Cols(child) ]),
	newType := TempArrayType(child, y, x, 1),
    StripList(
      List([ 1 .. Length(cols) ], (i) -> TempVec(TArray(TempArrayType(child,
               y, x, i), cols[i])))) );

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("~/Workspace/groupsSB/epi/A2/adModule.gi");
#

type:="A";
rank:=2;
nr_pos_roots:=3;
L:= SimpleLieAlgebra(type,rank,Rationals);
cb0:=CanonicalBasis(L);

M:=AdjointModule(L);
extM2:=ExteriorPowerOfAlgebraModule(M,2);

ade:=function(e)
	local mb,result,v;

    mb:=Basis(M);

	result:=[];
	for v in mb do
		Append(result,[Coefficients(mb,e^v)]);
	od;
	#result:=TransposedMat(result);
	return result;
end;

root_group:=function(index,t)
	local e,tmp,result,i;

    e:=ade(cb0[index]);
    tmp:=ShallowCopy(e);
	result:=tmp^0;

	i:=1;
	while Length(Set(Concatenation(tmp)))<>1 do
		result:=result+t^i*tmp/Factorial(i);
		i:=i+1;
		tmp:=tmp*e;
	od;
	return result;
end;

ad_nil_deg:=function(e)
    local ee,tmp;
    ee:=ade(e);
    tmp:=ShallowCopy(ee);
    i:=1;
    while Length(Set(Concatenation(tmp)))<>1 do;
        tmp:=tmp*ee;
        i:=i+1;
    od;
    return i;
end;

ad_exp:=function(e,t,v)
    local result,ord_e,i,j,tmp;

    ord_e:=ad_nil_deg(e);

    result:=v;
    for i in [1..ord_e-1] do
        tmp:=ShallowCopy(v);
        for j in [1..i] do
            tmp:=e^tmp;
        od;
        result:=result+t^i*tmp/Factorial(i);
    od;
    return result;
end;

diag_g:=function(e,wedge_sum)
    local fam_algmod,fam_wedge,o,result,new_wedge,i,bv,exp_term;
    fam_algmod:=FamilyObj(wedge_sum);
    o:=ExtRepOfObj(wedge_sum);
    fam_wedge:=FamilyObj(o);
    o:=ExtRepOfObj(o);

    result:=[];
    for i in [1,3..Length(o)-1] do
        new_wedge:=[];
        for bv in o[i] do
            exp_term:=ad_exp(e,1,bv);
            Add(new_wedge,exp_term);
        od;
        Add(result,new_wedge);
        Add(result,o[i+1]);
    od;

    result:=ObjByExtRep(fam_wedge,result);
    result:=ConvertToNormalFormMonomialElement(result);
    result:=ObjByExtRep(fam_algmod,result);

    return result;
end;

exp_action_mat_g:=function(e,M)
    local mb,v,result;
    mb:=Basis(M);

    result:=[];
 	for v in mb do
		Append(result,[Coefficients(mb,e^v)]);
	od;
	#result:=TransposedMat(result);
	return result;
end;

u1:=root_group(1,1);

f:=LeftModuleHomomorphismByMatrix(Basis(M),u1,Basis(M));
test:=Image(f,Basis(M)[1]);

bbb:=Basis(extM2);
tmp:=ExtRepOfObj(ExtRepOfObj(bbb[1]));
tmp:=ShallowCopy(tmp);
tmp:=[[Image(f,tmp[1][1]),Image(f,tmp[1][2])],1];

fam:=FamilyObj(ExtRepOfObj(bbb[1]));
tmp:=ObjByExtRep(fam,tmp);
tmp:=ConvertToNormalFormMonomialElement(tmp);

fam2:=FamilyObj(bbb[1]);
tmp:=ObjByExtRep(fam2,tmp);

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

#F ==========================================================================
#F VBlk(<M>, <v>)
#F
#F Matrix of <v> element diagonals -- each matrix entry is represented as a
#F list vector elements, or as vector values, e.g., TVect(t, v).value(..)
#F
Class(VBlk, SumsBase, Mat, rec(
    #-----------------------------------------------------------------------
    new := (self, M, v) >> SPL(WithBases(self, rec(
        element := M,
	v       := v,
        TType   := Cond( # NOTE: add checks to M
	    IsList(M),     TVect(UnifyTypes(List(Concatenation(M), InferType)).t, v),
	    IsSymbolic(M), M.t.t.t),
    ))).setDims(),
    #-----------------------------------------------------------------------
    dims := self >> Dimensions(self.element) * self.v,
    #-----------------------------------------------------------------------
    rChildren := self >> [self.element, self.v],
    rSetChild := rSetChildFields("element", "v"),
    #-----------------------------------------------------------------------
    toAMat   := self >> BlockMat(MapMat(self.element,
        x -> Diag(FList(self.TType.t, _unwrap(When(IsList(x), x, x.ev())))))).toAMat(),
    #-----------------------------------------------------------------------
    velement := self >> Cond(
	IsList(self.element), let(
	    d := Dimensions(self.element),
	    t := self.TType,
	    TArray(TArray(t, d[2]), d[1]).value(
		MapMat(self.element, x -> self.TType.value(x)))),
	IsSymbolic(self.element), self.element),
    #-----------------------------------------------------------------------
    transpose := self >> Cond(self.isSymmetric(), self,
        CopyFields(self, rec(element := TransposedMat(self.element),
                             dimensions := Reversed(self.dimensions)))),
    conjTranspose := self >> self.transpose(), # NOTE: implement this properly
    #-----------------------------------------------------------------------
    isSymmetric := False,
    setSymmetric := self >> Inherit(self, rec(isSymmetric := True)),
    #-----------------------------------------------------------------------
    isReal := self >> Cond(
	IsList(self.element),  ForAll(self.element, r -> ForAll(r, IsRealNumber)),
	not IsComplexT(self.element.t.t.t)),
));

#F ==========================================================================
#F RCVBlk(<M>, <v>)
#F
#F Matrix of complex diagonals --
#F     each matrix entry is represented as a list of v/2 complex numbers
#F
#F <v> here represents the *real* vector length
#F (YSV: it used to be complex vector length in older versions)
#F
Class(RCVBlk, VBlk, rec(
    #-----------------------------------------------------------------------
    toAMat := self >> BlockMat(MapMat(self.element,
        x -> RCDiag(RCData(FList(TComplex, x))))).toAMat(),
    #-----------------------------------------------------------------------
    velement := self >> Cond(
        IsList(self.element),
            apack.fromMat(self.element, cel -> ApplyFunc(vpack,
                ConcatList( Flat([cel]), e -> [re(e), im(e)] ))),
        IsSymbolic(self.element),
            self.element),
    #-----------------------------------------------------------------------
    isReal := self >> true,
    #-----------------------------------------------------------------------
    # transposition will conjugate complex numbers, similarly to RCDiag
    transpose := self >> ObjId(self)(
	MapMat(TransposedMat(self.element), x -> List(x, conj)), self.v),
    #-----------------------------------------------------------------------
    conjTranspose := self >> ObjId(self)(TransposedMat(self.element), self.v)
));

#F ==========================================================================
#F VReplicate(<v>) - 1 -> v vector,  each component is splatted.
#F
#F Takes 1 vector and produces <v> vectors, by vdup'ping (splatting) each component
#F
Class(VReplicate, BaseContainer, rec(
    #-----------------------------------------------------------------------
    new := (self, v) >> SPL(WithBases(self, rec(_children := [], v := v))).setDims(),
    #-----------------------------------------------------------------------
    dims := self >> [self.v^2, self.v],
    #-----------------------------------------------------------------------
    toAMat := self >> Tensor(
        I(self.v),
        Mat(TransposedMat([Replicate(self.v, 1)]))).toAMat(),
    #-----------------------------------------------------------------------
    print := (self,i,is) >> Print(self.name, "(", self.v, ")"),
    #-----------------------------------------------------------------------
    sums  := self >> self,
    #-----------------------------------------------------------------------
    isPermutation := False,
    #-----------------------------------------------------------------------
    rChildren := self >> [],
    rSetChild := (self, n, what) >> Error("no children"),
    #-----------------------------------------------------------------------
    transpose     := self >> InertTranspose(self),
    conjTranspose := self >> InertConjTranspose(self)
));

#F ==========================================================================
#F VHAdd(<v>) - takes a <v> vectors and produces 1 vector
#F                         - each vector is added horizontally
#F
Class(VHAdd, BaseContainer, rec(
    #-----------------------------------------------------------------------
    new := (self, v) >> SPL(WithBases(self, rec(v := v))).setDims(),
    #-----------------------------------------------------------------------
    dims := self >> [self.v, self.v^2],
    #-----------------------------------------------------------------------
    toAMat := self >> Tensor(
        I(self.v),
        Mat([Replicate(self.v,1)])).toAMat(),
    #-----------------------------------------------------------------------
    print := (self,i,is) >> Print(self.name, "(", self.v, ")"),
    #-----------------------------------------------------------------------
    sums := self >> self,
    #-----------------------------------------------------------------------
    isPermutation := False,
    #-----------------------------------------------------------------------
    rChildren := self >> [],
    _children := self >> [],
    #-----------------------------------------------------------------------
    transpose     := self >> Error("Not implemented"),
    conjTranspose := self >> Error("Not implemented"),
));


text:
--
-- In a parser using this template, the following macro may be redefined:
--
--     %additional_interfaces
--     %ast_class
--
-- B E G I N N I N G   O F   T E M P L A T E   btParserTemplateF
--
%Options programming_Language=typescript,margin=4,backtrack
%Options table,error_maps,scopes
%options prefix=TK_
%options action-block=("*.ts", "/.", "./")
%options ParseTable=ParseTable
%options nt-check

--
-- This template requires that the name of the EOF token be set
-- to EOF_TOKEN to be consistent with LexerTemplateD and LexerTemplateE
--
%EOF
    EOF_TOKEN
%End

%ERROR
    ERROR_TOKEN
%End

%Define

    $Header
    /.
                //
                // Rule %rule_number:  %rule_text
                //
                ./

    $BeginAction
    /.%Header%case %rule_number: {
                   //#line %next_line "%input_file%"./

    $EndAction
    /.            break;
                }./

    $BeginJava
    /.%Header%case %rule_number: {
                    %symbol_declarations
                    //#line %next_line "%input_file%"./

    $EndJava /.%EndAction./

    $NoAction
    /.%Header%case %rule_number:
                    break;./

    $BadAction
    /.%Header%case %rule_number:
                    throw ("No action specified for rule " + %rule_number);./

    $NullAction
    /.%Header%case %rule_number:
                    this.setResult(null);
                    break;./

    $BeginActions
    /.

        public  ruleAction(ruleNumber : number) : void
        {
            switch (ruleNumber)
            {./

    $SplitActions
    /.
                    default:
                        this.ruleAction%rule_number(ruleNumber);
                        break;
                }
                return;
            }

            public  ruleAction%rule_number(ruleNumber  : number) : void
            {
                switch (ruleNumber)
                {
                    //#line %next_line "%input_file%"./

    $EndActions
    /.
                default:
                    break;
            }
            return;
        }./

    $entry_declarations
    /.

        public  parse%entry_name(monitor? : Monitor, error_repair_count : number = 0) : %ast_class | null
        {
            this.btParser.setMonitor(monitor);

            try
            {
                return <%ast_class> this.btParser.fuzzyParseEntry(%sym_type.%entry_marker, error_repair_count);
            }
            catch (ex)
            {
                if( ex instanceof BadParseException ){
                    let e = <BadParseException>(ex);

                    this.prsStream.reset(e.error_token); // point to error token
                    let diagnoseParser = new DiagnoseParser(this.prsStream, %action_type.prsTable);
                    diagnoseParser.diagnoseEntry(%sym_type.%entry_marker, e.error_token);
                }
                else{
                    throw ex;
                }
            }

            return null;
        }
    ./

    --
    -- Macros that may be needed in a parser using this template
    --
    $additional_interfaces /../
    $ast_class /.%ast_type./
    $super_class /.Object./
    $unimplemented_symbols_warning /.false./

    --
    -- Old deprecated macros that should NEVER be used.
    --
    $setSym1 /. // macro setSym1 is deprecated. Use function setResult
                this.getParser().setSym1./
    $setResult /. // macro setResult is deprecated. Use function setResult
                 this.getParser().setSym1./
    $getSym /. // macro getSym is deprecated. Use function getRhsSym
              this.getParser().getSym./
    $getToken /. // macro getToken is deprecated. Use function getRhsTokenIndex
                this.getParser().getToken./
    $getIToken /. // macro getIToken is deprecated. Use function getRhsIToken
                 this.prsStream.getIToken./
    $getLeftSpan /. // macro getLeftSpan is deprecated. Use function getLeftSpan
                   this.getParser().getFirstToken./
    $getRightSpan /. // macro getRightSpan is deprecated. Use function getRightSpan
                    this.getParser().getLastToken./
%End

%Globals
    /.
import {BadParseException, RuleAction, PrsStream, ParseTable, BacktrackingParser, IToken, ErrorToken, ILexStream, NullExportedSymbolsException,
UnimplementedTerminalsException, Lpg, UndefinedEofSymbolException, NotBacktrackParseTableException, BadParseSymFileException,
IPrsStream, Monitor, DiagnoseParser, IAst, IAstVisitor, IAbstractArrayList, NotDeterministicParseTableException,
 DeterministicParser, NullTerminalSymbolsException } from "lpg2ts";
import { %prs_type } from ".\/%prs_type";
import { %sym_type } from ".\/%sym_type";
    ./
%End

%Headers
    /.
    export class %action_type extends %super_class implements RuleAction%additional_interfaces
    {
        private  prsStream  : PrsStream = new PrsStream();

        private  unimplementedSymbolsWarning : boolean = %unimplemented_symbols_warning;

        private static  prsTable : ParseTable = new %prs_type();
        public  getParseTable() : ParseTable { return %action_type.prsTable; }

        private  btParser : BacktrackingParser ;
        public  getParser() : BacktrackingParser{ return this.btParser; }

        private  setResult(object1 : any) : void{ this.btParser.setSym1(object1); }
        public  getRhsSym(i : number) : any{ return this.btParser.getSym(i); }

        public  getRhsTokenIndex(i : number) : number{ return this.btParser.getToken(i); }
        public  getRhsIToken(i : number) : IToken { return this.prsStream.getIToken(this.getRhsTokenIndex(i)); }

        public  getRhsFirstTokenIndex(i : number) : number { return this.btParser.getFirstToken(i); }
        public  getRhsFirstIToken(i : number) : IToken{ return this.prsStream.getIToken(this.getRhsFirstTokenIndex(i)); }

        public  getRhsLastTokenIndex(i : number):number { return this.btParser.getLastToken(i); }
        public  getRhsLastIToken(i : number):IToken { return this.prsStream.getIToken(this.getRhsLastTokenIndex(i)); }

        public getLeftSpan() :number { return this.btParser.getFirstToken(); }
        public  getLeftIToken() :IToken { return this.prsStream.getIToken(this.getLeftSpan()); }

        public getRightSpan() : number { return this.btParser.getLastToken(); }
        public  getRightIToken() : IToken { return this.prsStream.getIToken(this.getRightSpan()); }

        public  getRhsErrorTokenIndex(i : number) : number
        {
            let index = this.btParser.getToken(i);
            let err = this.prsStream.getIToken(index);
            return (err instanceof ErrorToken ? index : 0);
        }
        public  getRhsErrorIToken(i : number) : ErrorToken
        {
            let index = this.btParser.getToken(i);
            let err = this.prsStream.getIToken(index);
            return <ErrorToken> (err instanceof ErrorToken ? err : null);
        }

        public  reset(lexStream : ILexStream) : void
        {
            this.prsStream.resetLexStream(lexStream);
            this.btParser.reset(this.prsStream);

            try
            {
                this.prsStream.remapTerminalSymbols(this.orderedTerminalSymbols(), %action_type.prsTable.getEoftSymbol());
            }
            catch (e)
            {
                if( e instanceof NullExportedSymbolsException){

                }
                else if( e instanceof UnimplementedTerminalsException){
                    if (this.unimplementedSymbolsWarning) {
                        let unimplemented_symbols = e.getSymbols();
                        Lpg.Lang.System.Out.println("The Lexer will not scan the following token(s):");
                        for (let i : number = 0; i < unimplemented_symbols.size(); i++)
                        {
                            let id = <number>unimplemented_symbols.get(i);
                            Lpg.Lang.System.Out.println("    " + %sym_type.orderedTerminalSymbols[id]);
                        }
                        Lpg.Lang.System.Out.println();
                    }
                }
                else if( e instanceof UndefinedEofSymbolException){
                    throw  (new UndefinedEofSymbolException
                        ("The Lexer does not implement the Eof symbol " +
                        %sym_type.orderedTerminalSymbols[%action_type.prsTable.getEoftSymbol()]));
                }

            }
        }

        constructor(lexStream? :ILexStream)
        {
            super();
            try
            {
                this.btParser = new BacktrackingParser(null, %action_type.prsTable, <RuleAction> this);
            }
            catch (e)
            {
                if(e instanceof NotBacktrackParseTableException)
                throw (new NotBacktrackParseTableException
                                    ("Regenerate %prs_type.ts with -BACKTRACK option"));
                else if(e instanceof BadParseSymFileException){
                    throw (new BadParseSymFileException("Bad Parser Symbol File -- %sym_type.ts"));
                }
                else{
                    throw e;
                }
            }
            if(lexStream){
              this.reset(lexStream);
            }
        }

        public  numTokenKinds() :number { return %sym_type.numTokenKinds; }
        public  orderedTerminalSymbols()  : string[] { return %sym_type.orderedTerminalSymbols; }
        public  getTokenKindName(kind : number ) : string { return %sym_type.orderedTerminalSymbols[kind]; }
        public  getEOFTokenKind() : number{ return %action_type.prsTable.getEoftSymbol(); }
        public  getIPrsStream()  : IPrsStream{ return this.prsStream; }

        /**
         * @deprecated replaced by {@link #getIPrsStream()}
         *
         */
        public  getPrsStream()  : PrsStream{ return this.prsStream; }

        /**
         * @deprecated replaced by {@link #getIPrsStream()}
         *
         */
        public  getParseStream() : PrsStream { return this.prsStream; }

        public parser(error_repair_count : number = 0 ,  monitor? : Monitor) :  %ast_class | null
        {
            this.btParser.setMonitor(monitor);

            try
            {
                return <%ast_class> this.btParser.fuzzyParse(error_repair_count);
            }
            catch (ex)
            {
               if( ex instanceof BadParseException ){
                     let e = <BadParseException>(ex);
                    this.prsStream.reset(e.error_token); // point to error token

                    let diagnoseParser = new DiagnoseParser(this.prsStream, %action_type.prsTable);
                    diagnoseParser.diagnose(e.error_token);
                }
                else{
                    throw ex;
                }
            }

            return null;
        }

        //
        // Additional entry points, if any
        //
        %entry_declarations
    ./

%End

%Rules
    /.%BeginActions./
%End

%Trailers
    /.
        %EndActions
    }
    ./
%End

--
-- E N D   O F   T E M P L A T E
--


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

#F IsSPLMat(<x>) - returns true for parametrized matrices (Gath, Scat, Diag, etc)
#F
IsSPLMat := x -> IsRec(x) and IsBound(x._sym) and x._sym;

# ==========================================================================
# BaseMat
# ==========================================================================
Class(BaseMat, ClassSPL, rec(
    _mat := true,
    _short_print := true,
    #-----------------------------------------------------------------------
    isTerminal := True,
    isPermutation := False,
    children   := self >> [],
    terminate  := self >> self,
    # -------- Transformation rules support ---------------------------------
    rChildren := self >> [self.element],
    rSetChild := rSetChildFields("element"),
));


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

Declare(TestBench);

#F TestBench("name", <transforms> <opts>, <searchopts>)
#F
#F d := TestBench("default", [DFT(2), DFT(4)], SpiralDefaults, rec())
#F
#F TestBench Interface:
#F
#F     .build(transforms, opts, bopts, name, searchopts),
#F
#F     .generateCfiles := [true, false],  default = true
#F     .matrixVerify   := [true, false],  default = false
#F     .fftwVerify     := [true, false],  default = false
#F
#F     .fileTransform(t, opts),    .c filenames:   override in a subclass
#F     .funcTransform(t, opts),    function names: override in a subclass
#F     .txtFileName(runMethod),    timing file name
#F
#F     .run()              run DP
#F     .runExhaustive()    run exhaustive search
#F     .runRandom()        run 1 random tree, do not generate hash
#F     .runRandomSave()    run 1 random tree, generate and save hash
#F     .runRandom10()      run 10 random trees, do not generate hash
#F     .runRandomSave10()  run 10 random trees, generate and save hash
#F     .pickRandom()       generate random trees and save hash, do not run
#F
#F     .generateCode()
#F     .generateProductionCode()   same as .generateCode() but will use opts.production()
#F
#F     .entries()
#F     .times()
#F     .mflops()
#F
Class(TestBench, rec(
    ##
    ## Configuration
    ##
    generateCfiles    := true,
    matrixVerify      := false,
    fftwVerify        := false,
    outputExhaustive  := false,
    outputDir         := ".",
    fileTransform     := (self, t, opts) >> self.outputDir :: Conf("path_sep") :: self.name :: "_"
                                                   :: Drop(CodeletName(CodeletShape(t)), 1) :: ".c",
    funcTransform     := (self, t, opts) >> "sub1",
    prodFileTransform := (self, t, opts) >> self.fileTransform(t, opts), # used in generateProductionCode
    prodFuncTransform := (self, t, opts) >> self.funcTransform(t, opts),
    txtFileName       := (self, runMethod) >> self.outputDir :: Conf("path_sep") :: self.name :: "."
                                                  :: SubString(runMethod, 5) :: ".txt",
    ##
    ## Public methods
    ##

    # TestBench(<name>, <transforms>, <opts>, <search_opts>)
    __call__ := meth(self, name, transforms, opts, searchopts)
        local o;
	o := CopyFields(opts);
	if not IsBound(o.hashTable) then o.hashTable := HashTableDP(); fi;
	if not IsBound(o.hashFile)  then o.hashFile := Concat(self.outputDir, Conf("path_sep"), name, ".hash"); fi;
        return WithBases(self,
            rec(searchopts:=searchopts, opts:=o, transforms:=transforms, name:=name, verbosity:=1, callbacks:=[]));
    end,

    # TestBench.build(<transforms>, [<opts>], [<bench_opts>], [<bench_name>], [<search_opts>])
    #    Alternative constructor that supports default argument values if not provided
    #
    build := function(arg)
        local transforms, opts, bopts, name, dpr;

        transforms := When(IsList(arg[1]), arg[1], [arg[1]]);
        opts := When(Length(arg) >= 2, arg[2], SpiralDefaults);
        bopts := When(Length(arg) >= 3, arg[3], rec());
        name := When(Length(arg) >= 4, arg[4], "spiral");
        dpr := When(Length(arg) >= 5, arg[5], rec(verbosity := 0, timeBaseCases:=true));

        return CopyFields(TestBench(name, transforms, opts, dpr), bopts);
    end,

    generateCode           := self >> self._generateCode(self.transforms, self.opts),

    generateProductionCode := self >> self._generateCode(self.transforms, self.opts.production()),

    entries := self >> List(self.transforms, t -> self.entry(t)),

    entry := (self, t) >> let(l := self._rawentry(t),
	When(l=false, false, CopyFields(l, rec(ruletree := ApplyRuleTreeSPL(l.ruletree, t, self.opts))))),

    times   := self >> List(self.entries(), e -> e.measured),

    mflops  := self >> List(self.entries(), e -> self.acost(e) * LocalConfig.cpuinfo.freq / e.measured),

    # Use NonTerminal.normalizedArithCost() if it is there, otherwise return 0
    acost := (self, entry) >> self._rtflops(entry.ruletree),

    run            := arg >> arg[1]._run([],Drop(arg, 1), "_runDP",         true),
    runDP          := ~.run,
    runExhaustive  := arg >> arg[1]._run([], Drop(arg, 1), "_runExhaustive", true),
    runRandom      := arg >> arg[1]._run([arg[2]], Drop(arg, 2), "_runRandom",     false),
    runRandomSave  := arg >> arg[1]._run([arg[2]], Drop(arg, 2), "_runRandomSave", true),
    pickRandomSave := arg >> arg[1]._run([], Drop(arg, 1), "_pickRandomSave",true),

    ##
    ## Private methods
    ##

    _rtflops := rtree -> let(t:=rtree.node, When(IsBound(t.normalizedArithCost), EvalScalar(t.normalizedArithCost()), 0)),

    _verify := (self, opts, ruletree) >> VerifyMatrixRuleTree(ruletree, opts),

    # NOTE: Slightly hacked in. Check for opts.profile being bound etc. Look at VerifyMatrixRuleTree.
    _verifyfftw := (self, opts, code) >> opts.profile.verifyfftw(code, opts),

    _startHashFile := (self, hfile, d) >> PrintTo(hfile,
        "<# DPBench experiment '", self.name, "'\n",
        " # Started ", d[2], " ", d[3], " ", d[1], "  ", d[4], ":", d[5], "\n",
        " # Transforms: ", self.transforms, "#> \n\n",
        "ImportAll(spiral); Import(paradigms.common, paradigms.smp, platforms.sse, paradigms.vector); \n",
        "ImportAll(paradigms.vector); \n\n",
        "hash := HashTableDP(); \n"
    ),

    _loadHash := meth(self, hfile)
        local ns, result;
        ns := tab();
        result := READ(hfile, ns);
        if result = false or not IsBound(ns.hash) then return false;
        else return ns.hash;
        fi;
    end,

    _saveHash := meth(self, hfile, date, hash)
        local bucket, e;
        var.print := var.printFull;
        self._startHashFile(hfile, date);
        for bucket in hash.entries do
            for e in bucket do
                if e.data<>[] then
                    AppendTo(hfile, "HashAdd(hash, ", e.key, ", [", e.data[1], "]);\n");
                fi;
            od;
        od;
        var.print := var.printShort;
    end,

    reloadHash := meth(self)
       local hash, e;
       hash := self._loadHash(self.opts.hashFile);
       if (self.verbosity>0) then
	   PrintLine(When(hash=false, "Could not load ", "Loaded "), self.name, " (", self.opts.hashFile, ")");
       fi;
       if hash <> false then
	   self.opts.hashTable := hash;
       fi;
    end,

    _generateCode := meth(self, transforms, opts)
         local entries, e, c, t;
         for t in transforms do
             e := self.entry(t);
             if e = false then Error("Transform ", t, " not found in hashTable"); fi;
             c := CodeRuleTree(e.ruletree, opts);
             PrintLine(t, " -> ", self.prodFileTransform(t, opts));
             PrintTo(self.prodFileTransform(t, opts), PrintCode(self.prodFuncTransform(t, opts), c, opts));
         od;
    end,

    _showStats := meth(self, runMethod, t, ruletree, c, cycles, searchTime)
        local acc;
        # NOTE: slightly hacked in. Clean up to get both matrix and fftw verification to use already generated c.
	if self.matrixVerify or self.fftwVerify then
	    if self.matrixVerify then acc := self._verify(self.opts, ruletree);
	    else                      acc := self._verifyfftw(self.opts, c); fi;
	    _seqPerfStatsGflopsAcc(self.txtFileName(runMethod), t, self._rtflops(ruletree), cycles, searchTime, acc);
	else
	    When(self.opts.verbosity>-1,
		 _seqPerfStatsGflops(self.txtFileName(runMethod), t, LocalConfig.cpuinfo.freq, self._rtflops(ruletree), cycles, searchTime));
	fi;
    end,

    _rawentry := (self, t) >> let(
	lookup := MultiHashLookup(Concatenation([self.opts.hashTable], self.opts.baseHashes), HashAsSPL(t)),
	When(lookup=false or lookup=[], false, lookup[1])),

    #
    # Run methods
    #

    _runDP         := (self, t, opts) >> TimedAction(DP(t, self.searchopts, opts)),

    # Find best using an exhaustive search
    _runExhaustive := meth(self, t, opts)
       local r, searchTime, mincycles, mintree, rt, c, compiletime, cm, measuretime;

       r := AllRuleTrees(t, opts);
       searchTime := 0;
       mincycles := 10^100;

       for rt in r do
          [c,  compiletime] := TimedAction(CodeRuleTreeOpts(rt, opts));
          [cm, measuretime] := TimedAction(CMeasure(c, opts));
          if self.outputExhaustive then _seqPerfStatsGflops(
		  self.txtFileName("_runExhaustive-all"), t, LocalConfig.cpuinfo.freq, self._rtflops(rt), cm, compiletime+measuretime); fi;
          if (cm < mincycles) then
              mincycles := cm; mintree := Copy(rt);
          fi;
          searchTime:=searchTime+compiletime+measuretime;
       od;
       HashDelete(opts.hashTable,t);
       HashAdd(opts.hashTable, t, [rec(ruletree:=mintree, measured:=mincycles)]);

       return([mintree, searchTime]);
    end,

    # Run a <num> random ruletrees. Useful for quick, dirty, non-comprehensive tests.
    _runRandomNum := meth(self, num, t, opts)
       local r, c, start, cycles, i, res;
       start := TimeInSecs();
       res := [];
       for i in [1..num] do
           r := RandomRuleTree(t, opts);
           c := CodeRuleTree(r, opts);
           cycles := CMeasure(c, opts);
           Add(res, [r, cycles]);
       od;
       Sort(res, (a, b) -> a[2] < b[2]);
       return [res[1][1], res[1][2], TimeInSecs()-start];
     end,

    # Run random search and save result in hash.
    _runRandomSave := meth(self, num, t, opts)
       local r, srchTime, cycles;
       [r, cycles, srchTime] := self._runRandomNum(num, t, opts);
       HashDelete(opts.hashTable, t);
       HashAdd(opts.hashTable, t, [rec(ruletree:=r, measured:=cycles)]);
       return [r, srchTime];
    end,

    _runRandom := (self, num, t, opts) >> self._runRandomNum(num, t, opts){[1,3]},

    # Pick random and save result in hash (NOT measured).
    _pickRandomSave := meth(self, e, t, opts)
       local r, searchtime;
       [r, searchtime] := TimedAction(RandomRuleTree(t, opts));
       HashDelete(opts.hashTable,t);
       HashAdd(opts.hashTable,t,[rec(ruletree:=r)]);
       return [r, searchtime];
    end,

    _resume := self >> When(ForAny(self.entries(), e->e=false), self.reloadHash()),

    _run := meth(self, runArgs, transforms, runMethod, useHash)
        local t, outf, res, nopts, opts, c, cycles, hentry,  date, i, searchTime, acc, f, ruletree;
        MakeDir(self.outputDir);
        transforms := Flat(transforms);
        if transforms=[] then transforms := self.transforms; fi; # equivalent of runAll() in DPBench
        Constraint(ForAll(transforms, IsSPL));
        self._resume();
        for f in self.callbacks do f(self); od;
	opts := self.opts;
	date := Date();

	for t in transforms do
            # For run methods that use hash tables
	    if useHash then
		hentry := self._rawentry(t);
		if hentry = false then
		    res := ApplyFunc(self.(runMethod), runArgs :: [t, opts]);
		    if res[1] = [] then Error(runMethod, " did not find any ruletrees for <t> (", t, ")"); fi;
		    self._saveHash(opts.hashFile, date, opts.hashTable);
		    hentry := self._rawentry(t);
		    hentry.searchTime := res[2];
		    searchTime := res[2];
		else
		    searchTime := -1;
		fi;
		hentry.spectree := ApplyRuleTreeSPL(hentry.ruletree, t, opts);
                #NOTE: exhaustive search will not update cycles
		cycles := When(IsBound(hentry.measured), hentry.measured, 0);
		ruletree := hentry.spectree;
            # For run methods that don't use hash tables
	    else
		[ruletree, searchTime] := ApplyFunc(self.(runMethod), runArgs :: [t, opts]);
		cycles := 0;
	    fi;

            # HACK: It's a pain to do this in a cleaner way
	    if runMethod = "_pickRandomSave" then return; fi;

	    if self.generateCfiles then
                #NOTE: Shouldn't have to generate code or run this whole thing again.
		c := CodeRuleTree(ruletree, opts);
		cycles := CMeasure(c, opts);
		if useHash then hentry.measured := cycles; fi;

		nopts := CopyFields(opts, rec(fileinfo := rec(
			    cycles  := cycles,
			    flops   := self._rtflops(ruletree),
			    file    := self.fileTransform(t,opts),
			    algorithm := ruletree)));
		PrintTo(self.fileTransform(t,opts), PrintCode(self.funcTransform(t,opts), c, nopts));
            else
                c := skip();
	    fi;

	    self._showStats(runMethod, t, ruletree, c, cycles, searchTime);
        od;
    end,
));