text:
############################################################################
##
##  access.gd                     IRREDSOL                  Burkhard Höfling
##
##  Copyright © 2003–2016 Burkhard Höfling
##

############################################################################
##
#F  IndicesAbsolutelyIrreducibleSolubleMatrixGroups(<n>, <q>)
##
##  see the IRREDSOL manual
##
InstallGlobalFunction(IndicesAbsolutelyIrreducibleSolubleMatrixGroups,
    function(n, q)
        Info(InfoWarning, 1, "Obsolete function. See ?IndicesAbsolutelyIrreducibleSolubleMatrixGroups.");
        return IndicesIrreducibleSolubleMatrixGroups(n, q, 1);
    end);

############################################################################
##
#F  AbsolutelyIrreducibleSolubleMatrixGroup(<n>, <q>, <k>)
##
##  see the IRREDSOL manual
##
InstallGlobalFunction(AbsolutelyIrreducibleSolubleMatrixGroup,
    function(n, q, k)
        Info(InfoWarning, 1, "Obsolete function. See ? AbsolutelyIrreducibleSolubleMatrixGroup.");
        return IrreducibleSolubleMatrixGroup(n, q, 1, k);
    end);

############################################################################
##
#F  RecognitionAbsolutelyIrreducibleSolubleMatrixGroup(G, wantmat, wantgroup)
##
##  see the IRREDSOL manual
##
InstallGlobalFunction(RecognitionAbsolutelyIrreducibleSolubleMatrixGroup,
    function(G, wantmat, wantgroup)
        local r;
        Info(InfoWarning, 1, "Obsolete function. See ? RecognitionAbsolutelyIrreducibleSolubleMatrixGroup.");
        r := RecognitionIrreducibleSolubleMatrixGroup(G, wantmat, wantgroup);
        if r.id[3] <> 1 then
            Error("G is not absolutely irreducible");
        fi;
        r.id := r.id{[1,2,4]};
        return r;
    end);

############################################################################
##
#F  RecognitionAbsolutelyIrreducibleSolubleMatrixGroupNC(G, wantmat, wantgroup)
##
##  see the IRREDSOL manual
##
InstallGlobalFunction(RecognitionAbsolutelyIrreducibleSolubleMatrixGroupNC,
    function(G, wantmat, wantgroup)
        local r;
        Info(InfoWarning, 1, "Obsolete function. See ? RecognitionAbsolutelyIrreducibleSolubleMatrixGroupNC.");
        r := RecognitionIrreducibleSolubleMatrixGroupNC(G, wantmat, wantgroup);
        if r <> fail then
            if r.id[3] <> 1 then
                Error("G is not absolutely irreducible");
            fi;
            r.id := r.id{[1,2,4]};
        fi;
        return r;
    end);

############################################################################
##
#A  IdAbsolutelyIrreducibleSolubleMatrixGroup(<G>)
##
##  see the IRREDSOL manual
##
InstallGlobalFunction("IdAbsolutelyIrreducibleSolubleMatrixGroup",
    function(G)
        local r;
        Info(InfoWarning, 1, "Obsolete function. See ? IdAbsolutelyIrreducibleSolubleMatrixGroup.");
        r := IdIrreducibleSolubleMatrixGroup(G, false, false);
        if r[3] <> 1 then
            Error("G is not absolutely irreducible");
        fi;
        return r{[1,2,4]};
    end);

############################################################################
##
#E
##

max_stars_repo_name: fingolfin/irredsol

max_issues_repo_name: fingolfin/irredsol

max_forks_repo_name: fingolfin/irredsol


text:
--/**
-- * <copyright>
-- *
-- * Copyright (c) 2010, 2009 IBM Corporation and others.
-- * All rights reserved.   This program and the accompanying materials
-- * are made available under the terms of the Eclipse Public License v2.0
-- * which accompanies this distribution, and is available at
-- * http://www.eclipse.org/legal/epl-v20.html
-- *
-- * Contributors:
-- *   See (or edit) Notice Declaration below
-- *
-- * </copyright>
-- */
--
-- The Complete OCL Lexer
--

%options escape=$
%options la=2
%options fp=OCLLexer,prefix=Char_
%options single-productions
%options noserialize
%options package=org.eclipse.ocl.xtext.essentialocl.parser
%options template=../lpg/LexerTemplateF.gi
%options filter=OCLKWLexer.gi
%options export_terminals=("OCLParsersym.java", "TK_")
%options include_directory="../lpg"

%Import
	EssentialOCLLexer.gi
%End

%Define

	--
	-- Definition of macro used in the included file LexerBasicMap.g
	-- We redefine that one defined by EssentialOCLLexer
	--
	$kw_lexer_class /.OCLKWLexer./

%End

%Notice
	/./**
 * Complete OCL Lexer
 * <copyright>
 *
 * Copyright (c) 2010, 2009 IBM Corporation and others.
 * All rights reserved.   This program and the accompanying materials
 * are made available under the terms of the Eclipse Public License v2.0
 * which accompanies this distribution, and is available at
 * http://www.eclipse.org/legal/epl-v20.html
 *
 * Contributors:
 *   IBM - Initial API and implementation
 *   E.D.Willink - Bug 292112, 292594
 *   Adolfo Sanchez-Barbudo Herrera (Open Canarias) - LPG v 2.0.17 adoption (242153)
 *   Adolfo Sanchez-Barbudo Herrera (Open Canarias) - Introducing new LPG templates (299396)
 $copyright_contributions
 *******************************************************************************/
	./
%End

%Export
	AT
	CARET
	CARETCARET
	QUESTIONMARK
%End

%Rules
	Token ::= '@'
		/.$BeginAction
					makeToken($_AT);
		  $EndAction
		./

	Token ::= '^'
		/.$BeginAction
					makeToken($_CARET);
		  $EndAction
		./

	Token ::= '^' '^'
		/.$BeginAction
					makeToken($_CARETCARET);
		  $EndAction
		./

	Token ::= '?'
		/.$BeginAction
					makeToken($_QUESTIONMARK);
		  $EndAction
		./
%End

ext: gi


text:
# Calculate over the field Q(a) where a is a solution to x^2-3 = 0
x := Indeterminate(Rationals, "x");
p := x^2-3;
e := AlgebraicExtension(Rationals, p);
a := RootOfDefiningPolynomial(e);

b1 := [1, 0] * One(e);
b2 := [1/2, 1/2*a] * One(e);

s := [-1/2, 1/2*a] * One(e);
t := [-1/2, -1/2*a] * One(e);

u := s - b1;
v := t - b1;

ToFloat := function(element)
  local coefficients;
  coefficients := ExtRepOfObj(element);
  return Float(coefficients[1]) + Sqrt(3.0)*Float(coefficients[2]);
end;

scale := 4;
n := 3;
radius := 2*Ceil(n/2.0);

points := [];
for a in [-n..n] do
  for b in [-n..n] do
    for c in [b1, b2] do
      w := (a*u + b*v) + c;
      norm := ToFloat(w * w);
      if norm <= radius^2 then
        w := scale * w;
        Add(points, w);
      fi;
    od;
  od;
od;

PrintToConsole := function(index, element)
  Print(element, "\n");
end;

psFile := "data.ps";
PrintToPS := function(index, element)
  AppendTo(psFile, "\t[", ToFloat(element[1]), " ", ToFloat(element[2]), "]\n");
end;

tspFile := "data.tsp";
PrintToTSP := function(index, element)
  AppendTo(tspFile, index, " ", ToFloat(element[1]), " ", ToFloat(element[2]), "\n");
end;

header := Concatenation("NAME: hex-paper\n\
TYPE: TSP\n\
COMMENT: Taking concorde for a run on a small problem\n\
DIMENSION: ", String(Length(points)), "\n\
EDGE_WEIGHT_TYPE: EUC_2D\n\
NODE_COORD_TYPE: TWOD_COORDS\n\
NODE_COORD_SECTION:\n");

PrintTo(psFile);
PrintTo(tspFile, header);

index := 0;
for w in points do
  PrintToConsole(index, w);
  PrintToPS(index, w);
  PrintToTSP(index, w);
  index := index + 1;
od;

ext: gap

max_stars_repo_name: dvberkel/hexpaper-tsp

max_issues_repo_name: dvberkel/hexpaper-tsp

max_forks_repo_name: dvberkel/hexpaper-tsp


text:
# Built-in

[1, 3, -5]*[4, -2, -1];
# 3

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
############################################################################
##
##  primitive.gi                   IRREDSOL                 Burkhard Höfling
##
##  Copyright © 2003–2016 Burkhard Höfling
##

############################################################################
##
#F  PcGroupExtensionByMatrixAction(<pcgs>, <hom>)
##
##  Let <G> be a finite soluble group with pcgs <pcgs>, and let <hom> be a
##  group hom. $<hom>\colon G \to GL(n, p)$, where $p$ is a prime. Let  $E$
##  denote the split
##  extension of $G$ by $V = \F_p$, where <G> acts on <V> via <hom>.
##  This function returns a record with the following components.
##     ext:   the group $E$ as a new pc group
##     V:     the subgroup $V$ of $E$ corresponding to the vector space
##     C:     a complement of $V$ in $E$ isomorphic with $G$
##     embed: a group homomorphism $G \to E$ with image $C$
##     proj:  a group homomorphism $E \to G$ with kernel $V$
##     pcgsV: an induced pcgs of V (wrt. FamilyPcgs(E)) whose elements
##               correspond to the natural basis elements
##               of the vector space V
##     pcgsC: an induced pcgs of C (wrt. FamilyPcgs(E)) whose elements
##               correspond to the images of pcgs under embed;
##               the elements of pcgsC act on pcgsV as the images of pcgs
##               under hom act on the natural basis of V
##
InstallGlobalFunction(PcGroupExtensionByMatrixAction,
    function(pcgs, hom)
        local p, d, ros, f, coll, exp, mat, i, j, r, E, pcgsC, pcgsV;

        p := Size(FieldOfMatrixGroup(Range(hom)));
        d := DegreeOfMatrixGroup(Range(hom));
        if not IsPrimeInt(p) then
            Error("Range(hom) must be over a prime field ");
        fi;

        ros := RelativeOrders(pcgs);

        f := FreeGroup(Length(pcgs) + d);
        coll := SingleCollector(f,
            Concatenation(ros,
                ListWithIdenticalEntries(d, p)));

        # relations for complement - same as for those for G
        exp := [];
        exp{[1,3..2*Length(pcgs)-1]} := [1..Length(pcgs)];
        for i in [1..Length(pcgs)] do
            exp{[2,4..2*Length(pcgs)]} := ExponentsOfPcElement(pcgs, pcgs[i]^ros[i]);
            # Print("power relation ", i,": ", exp, "\n");
            SetPower(coll, i, ObjByExtRep(FamilyObj(f.1), exp));
            for j in [i+1..Length(pcgs)] do
                exp{[2,4..2*Length(pcgs)]} := ExponentsOfPcElement(pcgs, pcgs[j]^pcgs[i]);
                # Print("conj. relation ", j, "^", i,": ", exp, "\n");
            SetConjugate(coll, j, i, ObjByExtRep(FamilyObj(f.1), exp));
            od;
        od;

        # relations for socle
        for j in [1..d] do
            SetPower(coll, j+Length(pcgs), One(f));
        od;

        exp := [];
        exp{[1,3..2*d-1]} :=
            [Length(pcgs) + 1..Length(pcgs) + d];

        for i in [1..Length(pcgs)] do
            mat := ImageElm(hom, pcgs[i]);
            for j in [1..d] do
                exp{[2,4..2*d]} := List(mat[j], IntFFE);
                # Print("conj. relation ", j+ Length(pcgs), "^", i,": ", exp, "\n");
                SetConjugate(coll, j + Length(pcgs), i, ObjByExtRep(FamilyObj(f.1), exp));
            od;
        od;

        E := GroupByRwsNC(coll);
        SetSize(E, Product(ros) * p^d);
        pcgsV := InducedPcgsByPcSequenceNC(FamilyPcgs(E),
            FamilyPcgs(E){[Length(pcgs) + 1..Length(FamilyPcgs(E))]});

        # the following sets attributes/properties which are defined
        # in the CRISP packages

        pcgsC := InducedPcgsByPcSequenceNC(FamilyPcgs(E),
                FamilyPcgs(E){[1..Length(pcgs)]});

        r := rec(
            E := E,
            V := GroupOfPcgs(pcgsV),
            C := GroupOfPcgs(pcgsC),
            pcgsV := pcgsV,
            pcgsC := pcgsC);
        r.embed := GroupHomomorphismByImagesNC(GroupOfPcgs(pcgs), E, pcgs, pcgsC);
        SetIsInjective(r.embed, true);
        SetImagesSource(r.embed, r.C);
        r.proj := GroupHomomorphismByImagesNC(E, GroupOfPcgs(pcgs),
            Concatenation(pcgsC, pcgsV),
            Concatenation(pcgs, ListWithIdenticalEntries(d, OneOfPcgs(pcgs))));
        SetIsSurjective(r.proj, true);
        SetKernelOfMultiplicativeGeneralMapping(r.proj, r.V);
        return r;
    end);

############################################################################
##
#F  PrimitivePcGroupIrreducibleMatrixGroup(<G>)
##
##  see IRREDSOL documentation
##
InstallGlobalFunction(PrimitivePcGroupIrreducibleMatrixGroup,
    function(G)

        if not IsMatrixGroup(G) or not IsFinite(FieldOfMatrixGroup(G))
                or not IsPrimeInt(Size(FieldOfMatrixGroup(G)))
                or not IsIrreducibleMatrixGroup(G) then
            Error("G must be an irreducible matrix group over a prime field");
        fi;

        return PrimitivePcGroupIrreducibleMatrixGroupNC(G);
    end);

############################################################################
##
#F  PrimitivePcGroupIrreducibleMatrixGroupNC(<G>)
##
##  see IRREDSOL documentation
##
##  it is important that the map from Pcgs(Source(RepresentationIsomorphism(G)))
##
InstallGlobalFunction(PrimitivePcGroupIrreducibleMatrixGroupNC,
    function(G)

        local rep, ext;

        rep := RepresentationIsomorphism(G);
        ext := PcGroupExtensionByMatrixAction(Pcgs(Source(rep)), rep);
        SetSocle(ext.E, ext.V);
        SetSocleComplement(ext.E, ext.C);
        SetFittingSubgroup(ext.E, ext.V);

        # the following sets attributes/properties which are defined
        # in the CRISP packages

        if IsBoundGlobal("SetIsPrimitiveSoluble") then
            ValueGlobal("SetIsPrimitiveSoluble")(ext.E, true);
        fi;
        return ext.E;

    end);

############################################################################
##
#F  PrimitivePcGroup(<n>,<p>,<d>,<k>)
##
##  see IRREDSOL documentation
##
InstallGlobalFunction(PrimitivePcGroup,
    function(n, p, d, k)

        local q, desc, G, mat, bas, hom, ext, o, pcgs, pcgsC, pcgsV, H;

        if not IsPosInt(n) or not IsPosInt(d) or not IsPosInt(p) or not IsPrimeInt(p) or n mod d <> 0 then
            Error("n, p, and d must be positive integers, ",
                "p must be a prime, and d must divide n");
        elif not k in IndicesIrreducibleSolubleMatrixGroups(n, p, d) then
            Error("k must be in IndicesIrreducibleSolubleMatrixGroups(n, p, d)");
        else
            n := n /d;
            q := p^d;
            LoadAbsolutelyIrreducibleSolubleGroupData(n, q);
            if n > 1 then
                desc := IRREDSOL_DATA.GROUPS[n][q][k];
            fi;
            if not IsBound(IRREDSOL_DATA.PRIM_GUARDIANS[n]) then
                IRREDSOL_DATA.PRIM_GUARDIANS[n] := [];
            fi;
            if not IsBound(IRREDSOL_DATA.PRIM_GUARDIANS[n][q]) then
                if n = 1 then
                    G := IRREDSOL_DATA.GUARDIANS[1][q][1];
                    mat := [[Z(q)]];
                    hom := GroupHomomorphismByImagesNC(G, Group(mat),
                        MinimalGeneratingSet(G), [mat]);
                    SetIsBijective(hom, true);
                else
                    hom := IRREDSOL_DATA.GUARDIANS[n][q][desc[1]][3];
                    G := Source(hom);
                fi;
                if d > 1 then
                    bas := CanonicalBasis(AsVectorSpace(GF(p), GF(q)));
                    mat := List(InducedPcgsWrtFamilyPcgs(G),
                        g -> BlownUpMat(bas, ImageElm(hom, g)));
                    hom := GroupHomomorphismByImagesNC(G, Group(mat),
                        InducedPcgsWrtFamilyPcgs(G), mat);
                fi;

                IRREDSOL_DATA.PRIM_GUARDIANS[n][q] :=
                    PcGroupExtensionByMatrixAction(InducedPcgsWrtFamilyPcgs(G), hom);

            fi;
            ext := IRREDSOL_DATA.PRIM_GUARDIANS[n][q];
            if n = 1 then
                Assert(1, Length(MinimalGeneratingSet(ext.C)) = 1);
                o := IRREDSOL_DATA.GROUPS_DIM1[q][k][1];
                pcgsC := InducedPcgsByGenerators(FamilyPcgs(ext.E),
                    [MinimalGeneratingSet(ext.C)[1]^((q-1)/o)]);
            else
                pcgsC := CanonicalPcgsByNumber(ext.pcgsC, desc[2]);
            fi;
            pcgs := InducedPcgsByPcSequenceNC(FamilyPcgs(ext.E), Concatenation(pcgsC, ext.pcgsV));
            H := GroupOfPcgs(pcgs);
            SetIdPrimitiveSolubleGroup(H, [n*d,p,d,k]);
            SetSocle(H, ext.V);
            SetFittingSubgroup(H, ext.V);
            SetSocleComplement(H, GroupOfPcgs(pcgsC));

            # the following sets attributes/properties which are defined
            # in the CRISP packages

            if IsBoundGlobal("SetIsPrimitiveSoluble") then
                ValueGlobal("SetIsPrimitiveSoluble")(H, true);
            fi;
            return H;
        fi;
    end);

############################################################################
##
#F  IrreducibleMatrixGroupPrimitiveSolubleGroup(<G>)
##
##  see IRREDSOL documentation
##
InstallGlobalFunction(IrreducibleMatrixGroupPrimitiveSolubleGroup,
    function(G)

        local F, p, matgrp, compl;

        if not IsFinite(G) or not IsSolvableGroup(G) then
            Error("G must be finite and soluble");

        # test if primitive - use the CRISP method if it is available
        elif IsBoundGlobal("IsPrimitiveSoluble")
                 and ValueGlobal("IsPrimitiveSoluble")(G) then
            return IrreducibleMatrixGroupPrimitiveSolubleGroupNC(G);

        else # test for primitivity
            F := FittingSubgroup(G);

            if not IsPGroup(F)  or not IsAbelian(F) then
                Error("G must be primitive");
            else
                p := PrimePGroup(F);

                if ForAny(GeneratorsOfGroup(F), x -> x^p <> One(G)) then
                    Error("G must be primitive");
                else
                    matgrp := IrreducibleMatrixGroupPrimitiveSolubleGroupNC(G);
                    if not IsIrreducibleMatrixGroup(matgrp, GF(p)) then
                        Error("G must be primitive");
                    else
                        compl := ComplementClassesRepresentatives(G, F);
                        if Length(compl) <> 1 then
                            Error("G must be primitive");
                        fi;
                        SetSocle(G, F);
                        return matgrp;
                    fi;
                fi;
            fi;
        fi;
    end);

############################################################################
##
#F  IrreducibleMatrixGroupPrimitiveSolubleGroupNC(<G>)
##
##  see IRREDSOL documentation
##
InstallGlobalFunction(IrreducibleMatrixGroupPrimitiveSolubleGroupNC,
    function(G)

        local N, p, F, pcgsN, pcgsGmodN, GmodN, one, mat, mats, g, h, i, H, hom;

        N := FittingSubgroup(G);

        pcgsN := Pcgs(N);
        p := RelativeOrders(pcgsN)[1];
        F := GF(p);
        one := One(F);

        mats := [];

        pcgsGmodN := ModuloPcgs(G, N);
        for g in pcgsGmodN do
            mat := [];
            for i in [1..Length(pcgsN)] do
                mat[i] := ExponentsOfPcElement(pcgsN, pcgsN[i]^g)*one;
            od;
            Add(mats, ImmutableMatrix(F, mat));
        od;
        H := Group(mats);
        SetSize(H, Size(G)/Size(N));
        GmodN := PcGroupWithPcgs(pcgsGmodN);
        hom := GroupGeneralMappingByImagesNC(GmodN, H, FamilyPcgs(GmodN), mats);
        SetIsGroupHomomorphism(hom, true);
        SetIsBijective(hom, true);
        SetRepresentationIsomorphism(H, hom);
        return H;
    end);

############################################################################
##
#F  DoIteratorPrimitiveSolubleGroups(<convert_func>, <arg_list>)
##
##  generic constructor function for an iterator of all primitive soluble groups
##  which can construct permutation groups or pc groups (or other types of groups),
##  depending on convert_func
##
InstallGlobalFunction(DoIteratorPrimitiveSolubleGroups,
    function(convert_func, arg_list)

        local r, iter;

        r := CheckAndExtractArguments([
            [[Degree, NrMovedPoints, LargestMovedPoint], IsPosInt],
            [[Order, Size], IsPosInt]],
            arg_list,
            "IteratorPrimitivePcGroups");
        if ForAny(r.specialvalues, v -> IsEmpty(v)) then
            return Iterator([]);
        fi;

        iter := rec(convert_func := convert_func);

        if not IsBound(r.specialvalues[1]) then
            Error("IteratorPrimitivePcGroupsIterator: You must specify the degree(s) of the desired primitive groups");
        else
            iter.degs := Filtered(r.specialvalues[1], IsPPowerInt);
        fi;

        iter.degind := 0;

        if IsBound(r.specialvalues[2]) then
            iter.orders := r.specialvalues[2];
        else
            iter.orders := fail;
        fi;

        iter.iteratormatgrp := Iterator([]);

        iter.IsDoneIterator := function(iterator)

            local d, p, n, orders, o;

            if iterator!.degind > Length(iterator!.degs) then
                Error("isDoneIterator called after it returned true");
            fi;

            while IsDoneIterator(iterator!.iteratormatgrp) do
                iterator!.degind := iterator!.degind + 1;
                if iterator!.degind > Length(iterator!.degs) then
                    return true;
                fi;
                d := iterator!.degs[iterator!.degind];
                p := SmallestRootInt(d);
                n := LogInt(d, p);
                if IsAvailableIrreducibleSolubleGroupData(n, p) then
                    if iterator!.orders <> fail then
                        orders := [];
                        for o in iterator!.orders do
                            if o mod d = 0 then
                                Add(orders, o/d);
                            fi;
                        od;
                        iterator!.iteratormatgrp := IteratorIrreducibleSolubleMatrixGroups(
                            Degree, n, Field, GF(p), Order, orders);
                    else
                        iterator!.iteratormatgrp := IteratorIrreducibleSolubleMatrixGroups(
                            Degree, n, Field, GF(p));

                    fi;
                else
                    Error("groups of degree ", d, " are beyond the scope of the IRREDSOL library");
                    iterator!.iteratormatgrp := Iterator([]);
                fi;
            od;
            return false;
        end;

        iter.NextIterator := function(iterator)

            local G;

            G := NextIterator(iterator!.iteratormatgrp);
            return iterator!.convert_func(G);
        end;

        iter.ShallowCopy := function(iterator)
            return rec(
                orders := iterator!.orders,
                degs := iterator!.degs,
                degind := iterator!.degind,
                convert_func := iterator!.convert_func,
                iteratormatgrp := ShallowCopy(iterator!.iteratormatgrp),
                IsDoneIterator := iterator!.IsDoneIterator,
                NextIterator := iterator!.NextIterator,
                ShallowCopy := iterator!.ShallowCopy);
        end;
        return IteratorByFunctions(iter);
    end);

############################################################################
##
#F  IteratorPrimitivePcGroups(<func_1>, <val_1>, ...)
##
##  see the IRREDSOL manual
##
InstallGlobalFunction(IteratorPrimitivePcGroups,
    function(arg)
        return DoIteratorPrimitiveSolubleGroups(
            PrimitivePcGroupIrreducibleMatrixGroupNC,
            arg);
    end);

###########################################################################
##
#F  AllPrimitivePcGroups(<arg>)
##
##  see IRREDSOL documentation
##
InstallGlobalFunction(AllPrimitivePcGroups,
    function(arg)

        local iter, l, G;

        iter := CallFuncList(IteratorPrimitivePcGroups, arg);

        l := [];
        for G in iter do
            Add(l, G);
        od;
        return l;
    end);

###########################################################################
##
#F  OnePrimitivePcGroup(<arg>)
##
##  see IRREDSOL documentation
##
InstallGlobalFunction(OnePrimitivePcGroup,
    function(arg)

        local iter;

        iter := CallFuncList(IteratorPrimitivePcGroups, arg);
        if IsDoneIterator(iter) then
            return fail;
        else
            return NextIterator(iter);
        fi;
    end);

############################################################################
##
#F  PrimitivePermGroupIrreducibleMatrixGroup(<G>)
##
##  see IRREDSOL documentation
##
InstallGlobalFunction(PrimitivePermGroupIrreducibleMatrixGroup,
    function(G)

        if not IsMatrixGroup(G) or not IsFinite(FieldOfMatrixGroup(G))
                or not IsPrimeInt(Size(FieldOfMatrixGroup(G)))
                or not IsIrreducibleMatrixGroup(G) then
            Error("G must be an irreducible matrix group over a prime field");
        fi;

        return PrimitivePermGroupIrreducibleMatrixGroupNC(G);
    end);

############################################################################
##
#F  PrimitivePermGroupIrreducibleMatrixGroupNC(<G>)
##
##  see IRREDSOL documentation
##
InstallGlobalFunction(PrimitivePermGroupIrreducibleMatrixGroupNC,
    function( M )
        local  gensc, genss, V, bas, enum, G;
        V := FieldOfMatrixGroup( M ) ^ DimensionOfMatrixGroup( M );
        bas := CanonicalBasis(V);
        enum := EnumeratorByBasis(bas);
        gensc := List(GeneratorsOfGroup(M), x -> Permutation(x, enum));
        genss := List( bas, x -> Permutation( x, enum, \+));
        G := GroupByGenerators(Concatenation(genss, gensc));
        SetSize( G, Size( M ) * Size( V ) );
        SetSocle(G, Subgroup(G, genss));
        SetSocleComplement(G, Subgroup(G, gensc));

        # the following sets attributes/properties which are defined
        # in the CRISP packages

        if IsBoundGlobal("SetIsPrimitiveSoluble") then
            ValueGlobal("SetIsPrimitiveSoluble")(G, true);
        fi;
         return G;
    end);

############################################################################
##
#F  PrimitiveSolublePermGroup(<n>,<p>,<d>,<k>)
##
##  see IRREDSOL documentation
##
InstallGlobalFunction(PrimitiveSolublePermGroup,
    function(n, p, d, k)

        local G;
        if not IsPosInt(n) or not IsPosInt(d) or not IsPosInt(p) or not IsPrimeInt(p) then
            Error("n, p, and d must be positive integers, ",
                "p must be a prime, and d must divide n");
        elif not k in IndicesIrreducibleSolubleMatrixGroups(n, p, d) then
            Error("k must be in IndicesIrreducibleSolubleMatrixGroups(n, p, d)");
        else
            G := PrimitivePermGroupIrreducibleMatrixGroupNC(
                    IrreducibleSolubleMatrixGroup(n, p, d, k));
            SetIdPrimitiveSolubleGroup(G, [n,p,d,k]);
        fi;
        return G;
     end);

############################################################################
##
#F  IteratorPrimitiveSolublePermGroups(<func_1>, <val_1>, ...)
##
##  see the IRREDSOL manual
##
InstallGlobalFunction(IteratorPrimitiveSolublePermGroups,
    function(arg)
        return DoIteratorPrimitiveSolubleGroups(
            PrimitivePermGroupIrreducibleMatrixGroupNC,
            arg);
    end);

###########################################################################
##
#F  AllPrimitiveSolublePermGroups(<arg>)
##
##  see IRREDSOL documentation
##
InstallGlobalFunction(AllPrimitiveSolublePermGroups,
    function(arg)

        local iter, l, G;

        iter := CallFuncList(IteratorPrimitiveSolublePermGroups, arg);

        l := [];
        for G in iter do
            Add(l, G);
        od;
        return l;
    end);

###########################################################################
##
#F  OnePrimitiveSolublePermGroup(<arg>)
##
##  see IRREDSOL documentation
##
InstallGlobalFunction(OnePrimitiveSolublePermGroup,
    function(arg)

        local iter;

        iter := CallFuncList(IteratorPrimitiveSolublePermGroups, arg);
        if IsDoneIterator(iter) then
            return fail;
        else
            return NextIterator(iter);
        fi;
    end);

############################################################################
##
#E
##


text:
#
# cls 17: A_2
#
handle17char3:=function()
    local tmp1,tmp2,o,info,rep,cent,uu;
    o:=AllClasses(orbs)[17];
    info:=infos[17];

    rep:=Representative(o);
    #pr[15]=[0,0,1,1,1,0]
    #pr[13]=[0,1,1,1,0,0]
    #pr[26]=[1,1,1,2,1,0]
    tmp1:=ApplyRootsReflections(rep,[15,13,26]);
    uu:=Unipotent(chevalleyAdj(rep),[[3,1]]);
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[3],[-1]);

    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[15,13,26]);
    uu:=Unipotent(chevalleyAdj(cent),[[3,1]]);
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[3],[-1]);

    return [tmp1,tmp2];
end;

#
# cls 10: D_4(a_1)
#
handle10char3:=function()
    local tmp1,tmp2,tmp3,tmp4,o,info,rep,cent,uu;
    o:=AllClasses(orbs)[10];
    info:=infos[10];

    rep:=Representative(o);
    #pr[4] =[0,0,0,1,0,0]
    #pr[32]=[1,1,2,2,1,1]
    #pr[33]=[1,1,1,2,2,1]
    tmp1:=ApplyRootsReflections(rep,[4,32,33]);
    uu:=Unipotent(chevalleyAdj(rep),[[3,-1],[5,1]]);
    tmp1:=Conj(tmp1,uu);
    tmp1:=ConjugateByTori(tmp1,[4,1,2,6],[-1,-1,-1,-1]);

    #pr[20]=[0,1,0,1,1,1]
    #pr[21]=[0,0,1,1,1,1]
    tmp3:=ApplyRootsReflections(rep,[20,21]);
    uu:=Unipotent(chevalleyAdj(rep),[[3,1],[4,1],[5,-1],[9,1],[10,-1]]);
    tmp3:=Conj(tmp3,uu);
    tmp3:=ConjugateByTori(tmp3,[2,4],[-1,-1,-1,-1]);

    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ApplyRootsReflections(cent,[4,32,33]);
    uu:=Unipotent(chevalleyAdj(cent),[[3,-1],[5,1]]);
    tmp2:=Conj(tmp2,uu);
    tmp2:=ConjugateByTori(tmp2,[4,1,2,6],[-1,-1,-1,-1]);

    tmp4:=ApplyRootsReflections(cent,[20,21]);
    uu:=Unipotent(chevalleyAdj(cent),[[3,1],[4,1],[5,-1],[9,1],[10,-1]]);
    tmp4:=Conj(tmp4,uu);
    tmp4:=ConjugateByTori(tmp4,[2,4],[-1,-1,-1,-1]);

    return [tmp1,tmp3,tmp2,tmp4];
end;

#
# cls 4: E_6(a_3) [13,1,15,6,14,4]
#
handle4char3:=function()
    local tmp1,tmp2,o,info,rep,cent;
    o:=AllClasses(orbs)[4];
    info:=infos[4];

    rep:=Representative(o);
    tmp1:=ConjugateByTorus(rep,4,-1);# h_1 --> h_4

    cent:=FromPositiveBorel(o,info[2]);
    tmp2:=ConjugateByTorus(cent,4,-1);# h_1 --> h_4

    return [tmp1,tmp2];
end;


text:
merge_lists := function( list1, list2 )
  local i, j, merged_list;

  merged_list := [ ];

  i := 1;

  j := 1;

  while i <= Length( list1 ) and j <= Length( list2 ) do

      if list1[ i ] <= list2[ j ] then

          Add( merged_list, list1[ i ] );

          i := i + 1;

      else

          Add( merged_list, list2[ j ] );

          j := j + 1;

      fi;

  od;

  if i <= Length( list1 ) then

      Append( merged_list, list1{[ i .. Length( list1 ) ]} );

  else

      Append( merged_list, list2{[ j .. Length( list2 ) ]} );

  fi;

  return merged_list;

end;

merge_sort := function( list )
  local lists_to_merge;

  if Length( list ) = 1 then

      return list;

  fi;

  lists_to_merge := [ ];

  lists_to_merge[ 1 ] := merge_sort( list{[ 1 .. Int( Length( list ) / 2 ) ]} );

  lists_to_merge[ 2 ] := merge_sort( list{[ Int( Length( list ) / 2 ) + 1 .. Length( list ) ]} );

  return merge_lists( lists_to_merge[ 1 ], lists_to_merge[ 2 ] );

end;

max_stars_repo_name: Mynogs/Algorithm-Implementations

max_issues_repo_name: Mynogs/Algorithm-Implementations

max_forks_repo_name: Mynogs/Algorithm-Implementations


text:
############################################################################
##
##  recognizeprim.gi              IRREDSOL                  Burkhard Höfling
##
##  Copyright © 2003–2016 Burkhard Höfling
##

###########################################################################
##
#F  IdPrimitiveSolubleGroup(<grp>)
##
##  see IRREDSOL documentation
##
InstallMethod(IdPrimitiveSolubleGroup, "for soluble group",
     true, [IsSolvableGroup and IsFinite], 0,
     G -> IdIrreducibleSolubleMatrixGroup(IrreducibleMatrixGroupPrimitiveSolubleGroup(G)));

RedispatchOnCondition(IdPrimitiveSolubleGroup, true, [IsGroup],
     [IsFinite and IsSolvableGroup], 0);

###########################################################################
##
#F  IdPrimitiveSolubleGroupNC(<grp>)
##
##  see IRREDSOL documentation
##
InstallGlobalFunction(IdPrimitiveSolubleGroupNC,
     function(G)
          local id;
          id := IdIrreducibleSolubleMatrixGroup(IrreducibleMatrixGroupPrimitiveSolubleGroupNC(G));
          SetIdPrimitiveSolubleGroup(G, id);
          return id;
     end);

############################################################################
##
#F  RecognitionPrimitiveSolubleGroup(<G>)
##
##  see IRREDSOL documentation
##
InstallGlobalFunction(RecognitionPrimitiveSolubleGroup,
    function(G, wantiso)

        local N, F, p, pcgsN, C, pcgsC, one, i, mat, mats, CC, H, hom, infomat, info, rep, ext, imgs, g, r;

        N := FittingSubgroup(G);

        pcgsN := Pcgs(N);
        p := Set(RelativeOrders(pcgsN));
        if Length(p) <> 1 then
            Error("G must be primitive");
        fi;

        p := p[1];

        if not IsAbelian(N) or ForAny(pcgsN, g -> g^p <> One(G)) then
            Error("G must be primitive");
        fi;

        # now we know that N is elementary abelian of exponent p

        F := GF(p);
        one := One(F);

        mats := [];

        C := ComplementClassesRepresentatives(G, N);
        if Length(C) <> 1 then
          Error("G must be primitive");
        fi;

        C := C[1];

        # N is complemented

        pcgsC := Pcgs(C);
        for g in pcgsC do
            mat := [];
            for i in [1..Length(pcgsN)] do
                mat[i] := ExponentsOfPcElement(pcgsN, pcgsN[i]^g)*one;
            od;
            Add(mats, ImmutableMatrix(F, mat));
        od;

        if not MTX.IsIrreducible(GModuleByMats(mats, F)) then
            Error("G must be primitive");
        fi;

        H := Group(mats);

        # the recognition part works best if the source of the representation isomorphism is a pc group

        if IsPcGroup(C) then
            CC := C;
        else
            CC := PcGroupWithPcgs(Pcgs(C));
        fi;

        SetSize(H, Size(C));
        hom := GroupGeneralMappingByImagesNC(CC, H, Pcgs(CC), mats);
        SetIsGroupHomomorphism(hom, true);
        SetIsBijective(hom, true);
        SetRepresentationIsomorphism(H, hom);

        infomat := RecognitionIrreducibleSolubleMatrixGroup(H, wantiso, wantiso, wantiso);

        info := rec(id := infomat.id);
        if not wantiso then
            return info;
        fi;

        rep := RepresentationIsomorphism(infomat.group);

        ext := PcGroupExtensionByMatrixAction(Pcgs(Source(rep)), rep);

        imgs := [];
        for g in Pcgs(CC) do
            Add(imgs, ImageElm(ext.embed, ImageElm(infomat.iso, g)));
        od;
        for r in infomat.mat do
            g := PcElementByExponents(ext.pcgsV, List(r, IntFFE));
            Add(imgs, g);
        od;

        info.group := ext.E;
        info.iso := GroupHomomorphismByImages(G, ext.E,
            Concatenation(pcgsC, pcgsN), imgs);
        if info.iso = fail or not IsBijective(info.iso) then
            Error("wrong group homomorphism");
        fi;
        return info;
    end);

############################################################################
##
#E
##


text:
rot13 := function(s)
  local upper, lower, c, n, t;
  upper := "ABCDEFGHIJKLMNOPQRSTUVWXYZ";
  lower := "abcdefghijklmnopqrstuvwxyz";
  t := [ ];
  for c in s do
    n := Position(upper, c);
    if n <> fail then
      Add(t, upper[((n+12) mod 26) + 1]);
    else
      n := Position(lower, c);
      if n <> fail then
        Add(t, lower[((n+12) mod 26) + 1]);
      else
        Add(t, c);
      fi;
    fi;
  od;
  return t;
end;

a := "England expects that every man will do his duty";
# "England expects that every man will do his duty"
b := rot13(a);
# "Ratynaq rkcrpgf gung rirel zna jvyy qb uvf qhgl"
c := rot13(b);
# "England expects that every man will do his duty"


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

_dim3 := d -> "dim3("::d.id::")";

Class(HIPUnparser, CudaUnparser, rec(
    cu_call := (self, o, i, is) >>
        Print(Blanks(i), "hipLaunchKernelGGL(",
                self.infix([o.func, _dim3(o.dim_grid), _dim3(o.dim_block), "0", "0"]::o.args, ", "), ");\n"),
));

Class(FFTXHIPOpts, FFTXOpts, rec(
    tags := [],
    operations := rec(Print := s -> Print("<FFTX HIP options record>")),
    max_threads := 1024
));

doHIPify := function(opts)
    opts.originalCudaOptsID := Copy(opts.operations.Print);
    opts.operations := rec(Print := s -> Print("<FFTX HIPified CUDA options record>"));
    opts.unparser := HIPUnparser;
    opts.includes := ["\"hip/hip_runtime.h\""];
#    opts.postProcessCode := (c, opts) -> FixUpHIP_Code(PingPong_3Stages(c, opts), opts);
    opts.postProcessCode := (c, opts) -> FixUpHIP_Code(c, opts);

    return opts;
end;

Class(FFTXHIPDefaultConf, rec(
    __call__ := self >> self,
    getOpts := (self, t) >> doHIPify(ParseOptsCUDA(FFTXCUDADeviceDefaultConf, t)),
    operations := rec(Print := s -> Print("<FFTX FFTX HIPified CUDA Configuration>")),
    useHIP := true
));

hipConf := rec(
    defaultName := "defaultHIPConf",
    defaultOpts := (arg) >> FFTXHIPDefaultConf,
    confHandler := doHIPify
);

fftx.FFTXGlobals.registerConf(hipConf);

max_stars_repo_name: mikefranusich/spiral-package-fftx

max_issues_repo_name: mikefranusich/spiral-package-fftx

max_forks_repo_name: mikefranusich/spiral-package-fftx


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

# DSCirculant(<n>, <P>, <v>, <ds>, <ofs>)
#    Downsampled circulant
#
Class(DSCirculant, Circulant, rec(
    abbrevs := [ (n, P, v, ds, ofs) ->
	    [ Checked(IsPosInt(n), n),
	      toFunc(P),
	      Checked(IsInt(v), v),
	      Checked(IsPosInt(ds), ds),
	      Checked(IsPosInt0(ofs), ofs < ds, ofs) ]
    ],

    dims := self >> let(p:=self.params, n:=p[1], ds:=p[4], ofs:=p[5],
	[ FloorRat((n+ofs)/ds), n ]),

    terminate := self >> let(cterminate := Circulant.terminate,
        DownSample(self.params[1], self.params[4], self.params[5]) *
        cterminate(self)),

    transpose := NonTerminal.transpose,

    HashId := self >> let(p := self.params, [p[4], p[1], p[2].domain(), p[3]])
));

# since Circulant is already added, DSCirculant won't be added automatically
AddNonTerminal(DSCirculant);

##
## Rules
##
RulesFor(DSCirculant, rec(
    ###################################################################
    ## Time domain methods
    ###################################################################

    #F RuleCirculant_Base: (base case)
    #F
    DSCirculant_Base := rec (
	info             := "DSCirculant -> Mat",
	forTransposition := false,
	isApplicable     := P -> P[1] <= d_sqrt(SPL_DEFAULTS.globalUnrolling),
	allChildren      := P -> [[ ]],
	rule := (P, C) -> ApplyFunc(DSCirculant, P).terminate()
    ),

    DSCirculant_toFilt := rec (
	info             := "DSCirculant -> VStack(.., Filt, ...)",
	forTransposition := false,
	isApplicable     := P -> P[1] > 2 and (P[1]-P[2].domain()) > 2 and P[3] <= 0,
	allChildren      := P -> [[ ]],

	rule := function(P, C)
	    local n,l,r,coeffs,nc,ds,ofs,ofs_l,ofs_r,ofs_filt,ds_l,ds_r,ds_filt,j,jj,k,kk,v,
	         filtpoly,filtcoef,filtlen,f;

	    n := P[1];  coeffs := P[2]; nc := coeffs.domain();
	    l := -P[3]; r := -l + nc - 1;

	    ds := P[4]; ofs := P[5];

	    ofs_l := ofs;
	    ofs_filt := (l - ofs_l) mod ds;
	    ofs_r := (n-l-r - ofs_filt) mod ds;

	    ds_l    := CeilingRat((l - ofs_l) / ds);
	    ds_filt := CeilingRat((n-l-r - ofs_filt) / ds);
	    ds_r    := CeilingRat((r - ofs_r) / ds);

	    j := Ind( ds_l );
	    k := Ind( ds_r );
	    f := Ind( ds_filt );
	    v := Ind(nc);
	    jj := Ind(l);
	    kk := Ind(r);

	    filtpoly := Poly(List(coeffs.tolist(),EvalScalar), -l);
	    filtcoef := FillZeros(filtpoly)[1];
	    filtlen := Length(filtcoef);

	    return
	    VStack(
	        BB(ISum(j, j.range,
		     Data(jj, ds*j+ofs_l,
			 Scat(fBase(j)) *
			 RowVec(fCompose(coeffs, Lambda(v, cond(leq(v,r+jj), v-jj+l, v-jj+l-nc)))) *
			 Gath(n, nc, Lambda(v, cond(leq(v,r+jj), v, v+(n-l-r)-1)))))),

		ISum(f, f.range,
		    Scat(fBase(f)) *
		    Blk([ filtcoef ]) *
		    Gath(fAdd(n, filtlen, ds*f + ofs_filt))),

		BB(ISum(k, k.range,
		     Data(kk, ds*k+ofs_r,
			 Scat(fBase(k)) *
			 RowVec(fCompose(coeffs, Lambda(v, cond(leq(v,kk), v-kk+(nc-1), v-kk-1)))) *
			 Gath(n, nc, Lambda(v, cond(leq(v,kk), v, v+(n-l-r)-1))))))
	    );
	end
    )
));


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

_dim3 := d -> "dim3("::d.id::")";

Class(HIPUnparser, CudaUnparser, rec(
    cu_call := (self, o, i, is) >>
        Print(Blanks(i), "hipLaunchKernelGGL(",
                self.infix([o.func, _dim3(o.dim_grid), _dim3(o.dim_block), "0", "0"]::o.args, ", "), ");\n"),
));

max_stars_repo_name: franzfranchetti/spiral-package-fftx

max_issues_repo_name: franzfranchetti/spiral-package-fftx

max_forks_repo_name: franzfranchetti/spiral-package-fftx


text:
# Using SageMath interface for GAP
# To start GAP inside of SageMath:
# gap.console()
#
# Defining groups in GAP:
# PSL := FreeProduct(CyclicGroup(2),CyclicGroup(3))
#
# Using SageMath to access GAP's
# finitely presented groups functionality
# https://doc.sagemath.org/html/en/reference/groups/sage/groups/finitely_presented.html

# from sage.interfaces.gap import get_gap_memory_pool_size, set_gap_memory_pool_size
# set_gap_memory_pool_size(20000000000)

k := 4;

F := FreeGroup("t_k", "c", "d");

G := F / [F.3^(-1)*F.1*(F.3*F.2^(-1))^(-k+1)*F.1*(F.3*F.2^(-1))^(k-1)*(F.1^(-1)),
          F.2^(-1)*F.1*F.2*F.3^(-1)*(F.1)^(-1)*F.3*F.1*F.3*F.2^(-1)*(F.1)^(-1)];

RelatorsOfFpGroup(G);

AllHomomorphismClasses(G, SymmetricGroup(3));
# AllHomomorphisms(G, SymmetricGroup(3));


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

Declare(ISRDFT, ISRDFT2, ISRDFT3, ISRDFT4);

Class(SRDFT, PRDFT, rec(
    dims := self >> [self.params[1], self.params[1]],
    inverse := self >> ISRDFT(self.params[1], -self.params[2]).setTransposed(self.transposed),
    inverseViaDiag := self >> self.transpose() * ISRDFT.diag(self.params[1]),
    terminate := self >> let(N := self.params[1], k := self.params[2],
	m  := MatSPL(PRDFT(N,k)),
	sm := Mat(Cond(IsOddInt(N),
	         Concatenation([m[1]], m{[3..Length(m)]}),
	         Concatenation([m[1]], m{[3..Length(m)-2]}, [m[Length(m)-1]]))),
	Cond(self.transposed, sm.transpose(), sm)),
));
SRDFT1 := SRDFT;
Class(SRDFT2, PRDFT2, rec(
    dims := self >> [self.params[1], self.params[1]],
    inverse := self >> ISRDFT3(self.params[1], -self.params[2]).setTransposed(self.transposed),
    terminate := self >> let(N := self.params[1], k := self.params[2],
	m  := MatSPL(PRDFT2(N,k)),
	sm := Mat(Cond(IsOddInt(N),
	         Concatenation([m[1]], m{[3..Length(m)]}),
	         Concatenation([m[1]], m{[3..Length(m)-2]}, [m[Length(m)]]))),
	Cond(self.transposed, sm.transpose(), sm)),
));
Class(SRDFT3, PRDFT3, rec(
    dims := self >> [self.params[1], self.params[1]],
    inverse := self >> ISRDFT2(self.params[1], -self.params[2]).setTransposed(self.transposed),
    terminate := self >> let(N := self.params[1], k := self.params[2],
	m  := MatSPL(PRDFT3(N,k)),
	sm := Mat(Cond(IsOddInt(N),
	         m{[1..Length(m)-1]},
	         m)),
	Cond(self.transposed, sm.transpose(), sm)),
));
Class(SRDFT4, PRDFT4, rec(
    dims := self >> [self.params[1], self.params[1]],
    inverse := self >> ISRDFT4(self.params[1], -self.params[2]).setTransposed(self.transposed),
    terminate := self >> let(N := self.params[1], k := self.params[2],
	m  := MatSPL(PRDFT4(N,k)),
	sm := Mat(Cond(IsOddInt(N),
	         Concatenation(m{[1..Length(m)-2]}, [m[Length(m)]]),
	         m)),
	Cond(self.transposed, sm.transpose(), sm)),
));

Class(ISRDFT, IPRDFT, rec(
    dims := self >> [self.params[1], self.params[1]],
    diag := n -> When(n=2, I(2),
        Diag(diagDirsum(fConst(TReal, 1, 1), fConst(TReal, 2*Int((n-1)/2), 2),
                When(IsEvenInt(n), fConst(TReal, 1, 1), [])))),

    terminate := self >> let(N := self.params[1], k := self.params[2],
	m  := MatSPL(IPRDFT(N,k)),
	sm := Mat(Cond(IsOddInt(N),
	         List(m, r->Concatenation([r[1]], r{[3..Length(r)]})),
	         List(m, r->Concatenation([r[1]], r{[3..Length(r)-2]}, [r[Length(r)-1]])))),
	Cond(self.transposed, sm.transpose(), sm)),
));
Class(ISRDFT2, IPRDFT2, rec(
    dims := self >> [self.params[1], self.params[1]],
    terminate := self >> let(N := self.params[1], k := self.params[2],
	m  := MatSPL(IPRDFT2(N,k)),
	sm := Mat(Cond(IsOddInt(N),
	         List(m, r->r{[1..Length(r)-1]}),
	         m)),
	Cond(self.transposed, sm.transpose(), sm)),
));
Class(ISRDFT3, IPRDFT3, rec(
    dims := self >> [self.params[1], self.params[1]],
    terminate := self >> let(N := self.params[1], k := self.params[2],
	m  := MatSPL(IPRDFT3(N,k)),
	sm := Mat(Cond(IsOddInt(N),
	         List(m, r->Concatenation([r[1]], r{[3..Length(r)]})),
	         List(m, r->Concatenation([r[1]], r{[3..Length(r)-2]}, [r[Length(r)]])))),
	Cond(self.transposed, sm.transpose(), sm)),
));
Class(ISRDFT4, IPRDFT4, rec(
    dims := self >> [self.params[1], self.params[1]],
    terminate := self >> let(N := self.params[1], k := self.params[2],
	m  := MatSPL(IPRDFT4(N,k)),
	sm := Mat(Cond(IsOddInt(N),
	         List(m, r->Concatenation(r{[1..Length(r)-2]}, [r[Length(r)]])),
	         m)),
	Cond(self.transposed, sm.transpose(), sm)),
));

RulesFor(SRDFT1, rec(
    SRDFT1_toPRDFT1 := rec(
	isApplicable := P -> true,
	allChildren := P -> [[ PRDFT1(P[1], P[2]) ]],
	rule := (P, C) -> Pack_CCS(P[1]) * C[1]
    )
));

RulesFor(SRDFT3, rec(
    SRDFT3_toPRDFT3 := rec(
	isApplicable := P -> true,
	allChildren := P -> [[ PRDFT3(P[1], P[2]) ]],
	rule := (P, C) -> Pack_CCS3(P[1]) * C[1]
    )
));


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

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

################################
#         MACRO_2xf
#

Class(vunpacklo_2xf, _macro_binop_ev, VecExp_2.binary(), rec(
    semantic := (in1, in2, p) -> unpacklo(in1, in2, 2, 1)));

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

Class(vpacklo_2xf, _macro_binop_ev, VecExp_2.binary(), rec(
    semantic := (in1, in2, p) -> shuffle(in1, in2, [1, 1], 2, 1)));

Class(vpackhi_2xf, _macro_binop_ev, VecExp_2.binary(), rec(
    semantic := (in1, in2, p) -> shuffle(in1, in2, [2, 2], 2, 1)));

Class(vmulcx_2xf, VecExp_2.binary());
Class(vswapcx_2xf, VecExp_2.unary());

################################
#         MACRO_4xf
#

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

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

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

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

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

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

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

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

Class(vmulcx_4xf, VecExp_4.binary());
Class(vswapcx_4xf, VecExp_4.unary());

################################
#         MACRO_8xf
#

Class(vunpacklo_8xf, _macro_binop_ev, VecExp_8.binary(), rec(
    semantic := (in1, in2, p) -> unpacklo(in1, in2, 8, 1)));

Class(vunpackhi_8xf, _macro_binop_ev, VecExp_8.binary(), rec(
    semantic := (in1, in2, p) -> unpackhi(in1, in2, 8, 1)));

Class(vunpacklo2_8xf, _macro_binop_ev, VecExp_8.binary(), rec(
    semantic := (in1, in2, p) -> unpacklo(in1, in2, 8, 2)));

Class(vunpackhi2_8xf, _macro_binop_ev, VecExp_8.binary(), rec(
    semantic := (in1, in2, p) -> unpackhi(in1, in2, 8, 2)));

Class(vpacklo_8xf, _macro_binop_ev, VecExp_8.binary(), rec(
    semantic := (in1, in2, p) -> shuffle(in1, in2, [1, 3, 5, 7, 1, 3, 5, 7], 8, 1)));

Class(vpackhi_8xf, _macro_binop_ev, VecExp_8.binary(), rec(
    semantic := (in1, in2, p) -> shuffle(in1, in2, [2, 4, 6, 8, 2, 4, 6, 8], 8, 1)));

Class(vpacklo2_8xf, _macro_binop_ev, VecExp_8.binary(), rec(
    semantic := (in1, in2, p) -> shuffle(in1, in2, [1, 3, 1, 3], 8, 2)));

Class(vpackhi2_8xf, _macro_binop_ev, VecExp_8.binary(), rec(
    semantic := (in1, in2, p) -> shuffle(in1, in2, [2, 4, 2, 4], 8, 2)));

Class(vtrcx_8xf, VecExp_8.unary(), rec(
    semantic := (in1, p) -> shuffle(in1, in1, [1, 2, 5, 6, 3, 4, 7, 8], 8, 1),
    ev := self >> self.t.value(self.semantic(self._vval(1), [])).ev()
));

Class(vmulcx_8xf, VecExp_8.binary());
Class(vswapcx_8xf, VecExp_8.unary());

#################################
# Three Op version
#
ThreeOpFromBinOp(
    Class(vunpacklo_2xf_cmd,  assign_cmd, rec( exp_op := vunpacklo_2xf )),
    Class(vunpackhi_2xf_cmd,  assign_cmd, rec( exp_op := vunpackhi_2xf )),
    Class(vpacklo_2xf_cmd,    assign_cmd, rec( exp_op := vpacklo_2xf   )),
    Class(vpackhi_2xf_cmd,    assign_cmd, rec( exp_op := vpackhi_2xf   )),
    Class(vmulcx_2xf_cmd,     assign_cmd, rec( exp_op := vmulcx_2xf    )),

    Class(vunpacklo_4xf_cmd,  assign_cmd, rec( exp_op := vunpacklo_4xf  )),
    Class(vunpackhi_4xf_cmd,  assign_cmd, rec( exp_op := vunpackhi_4xf  )),
    Class(vunpacklo2_4xf_cmd, assign_cmd, rec( exp_op := vunpacklo2_4xf )),
    Class(vunpackhi2_4xf_cmd, assign_cmd, rec( exp_op := vunpackhi2_4xf )),
    Class(vpacklo_4xf_cmd,    assign_cmd, rec( exp_op := vpacklo_4xf    )),
    Class(vpackhi_4xf_cmd,    assign_cmd, rec( exp_op := vpackhi_4xf    )),
    Class(vpacklo2_4xf_cmd,   assign_cmd, rec( exp_op := vpacklo2_4xf   )),
    Class(vpackhi2_4xf_cmd,   assign_cmd, rec( exp_op := vpackhi2_4xf   )),
    Class(vmulcx_4xf_cmd,     assign_cmd, rec( exp_op := vmulcx_4xf     )),

    Class(vunpacklo_8xf_cmd,  assign_cmd, rec( exp_op := vunpacklo_8xf  )),
    Class(vunpackhi_8xf_cmd,  assign_cmd, rec( exp_op := vunpackhi_8xf  )),
    Class(vunpacklo2_8xf_cmd, assign_cmd, rec( exp_op := vunpacklo2_8xf )),
    Class(vunpackhi2_8xf_cmd, assign_cmd, rec( exp_op := vunpackhi2_8xf )),
    Class(vpacklo_8xf_cmd,    assign_cmd, rec( exp_op := vpacklo_8xf    )),
    Class(vpackhi_8xf_cmd,    assign_cmd, rec( exp_op := vpackhi_8xf    )),
    Class(vpacklo2_8xf_cmd,   assign_cmd, rec( exp_op := vpacklo2_8xf   )),
    Class(vpackhi2_8xf_cmd,   assign_cmd, rec( exp_op := vpackhi2_8xf   )),
    Class(vmulcx_8xf_cmd,     assign_cmd, rec( exp_op := vmulcx_8xf     ))
);

Class(vswapcx_4xf_cmd,    assign_cmd, rec( exp_op := vswapcx_4xf    ));
Class(vswapcx_2xf_cmd,    assign_cmd, rec( exp_op := vswapcx_2xf    ));
Class(vswapcx_8xf_cmd,    assign_cmd, rec( exp_op := vswapcx_8xf    ));

Class(vtrcx_8xf_cmd, assign_cmd, rec(exp_op := vtrcx_8xf));

Class(vstore_cmd, assign_cmd);
Class(vload_cmd,  assign_cmd);

RewriteRules(ThreeOpRuleSet, rec(
    assign_to_store := Rule([assign, @(1, nth), @(2)], e -> vstore_cmd(@(1).val.loc, @(1).val.idx, @(2).val)),
    assign_to_load  := Rule([assign, @(1, var), @(2, nth)], e -> vload_cmd(@(1).val, @(2).val.loc, @(2).val.idx)),
));


text:
# Apparently GAP can only remove a file, not a directory
RemoveFile("input.txt");
# true
RemoveFile("docs");
# fail


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

# from FFTX C++ unparser
Class(TFCall, Tagged_tSPL_Container, rec(
    abbrevs :=  [ (nt, cconf) -> Checked(IsSPL(nt), [nt, cconf]) ],
    transpose := self >> ObjId(self)(self.params[1].transpose(), self.params[2]).withTags(self.getTags())
));

Class(TDeviceCall, Tagged_tSPL_Container, rec(
    abbrevs :=  [ (nt, cconf) -> Checked(IsSPL(nt), [nt, cconf]) ],
    transpose := self >> ObjId(self)(self.params[1].transpose(), self.params[2]).withTags(self.getTags())
));

Class(TDecl, Tagged_tSPL_Container, rec(
    abbrevs :=  [ (nt, vlist) -> Checked(IsSPL(nt), [nt, vlist]) ],
    transpose := self >> ObjId(self)(self.params[1].transpose(), self.params[2]).withTags(self.getTags())
));

Class(TDAGNode, Tagged_tSPL_Container, rec(
    abbrevs :=  [ (nt, ylist, xlist) -> Checked(IsSPL(nt), [nt, When(IsList(ylist), ylist, [ylist]), When(IsList(xlist), xlist, [xlist])]) ],
    transpose := self >> ObjId(self)(self.params[1].transpose(), self.params[2], self.params[3]).withTags(self.getTags())
));

Class(TDAG, TCompose, rec(
    terminate := self >> Error("Not yet implemented."),

    from_rChildren := (self, rch) >> let(
        len := Length(rch),
        transposed := rch[len-1],
        tags := rch[len],
        t := ApplyFunc(ObjId(self), [rch{[1..len-2]}]),
        tt := When(transposed, t.transpose(), t),
        attrTakeA(tt.withTags(tags), self)
    ),

    rChildren := self >>
        Concatenation(self.params[1], [self.transposed, self.tags]),

    rSetChild := meth(self, n, newChild)
        local l;
        l := Length(self.params[1]);
        if n <= l then
            self.params[1][n] := newChild;
        elif n = l+1 then
            self.transposed := newChild;
        elif n = l+2 then
            self.tags := newChild;
        else Error("<n> must be in [1..", l+2, "]");
        fi;
        # self.canonizeParams(); ??
        self.dimensions := self.dims();
    end,

));

## the dims here are ordered as in the MDDFT case: e.g., for 3D: odims[1] = x, odims[2] = y, odims[3] = z
#Class(TResample, Tagged_tSPL_Container, rec(
#    abbrevs :=  [ (odims, idims) -> [ odims, idims, [] ],
#                  (odims, idims, shifts) -> [ odims, idims, shifts ],],
#    dims := self >> [Product(self.params[1]), Product(self.params[2])],
#    isReal := True,
#    terminate := self >> Error("Not yet implemented."),
#    transpose := self >> ObjId(self)(self.params[2], self.params[1], self.params[3]).withTags(self.getTags())
#));

Declare(TSparseMat);
Class(TSparseMat, Tagged_tSPL_Container, rec(
    abbrevs :=  [ (dims, entries) -> [ dims, entries ] ],
    dims := self >> self.params[1],
    terminate := self >> Error("Not yet implemented."),
    transpose := self >> TSparseMat(Reversed(self.params[1]), [x->Error("Not yet done")])
));

Declare(TIterVStack);

Class(TIterHStack, Tagged_tSPL_Container, rec(
    abbrevs :=  [ (nt, idx) -> Checked(IsSPL(nt), IsVar(idx),
	    [nt, idx]) ],

    dims := self >> self.params[1].dims(){[1]} :: self.params[1].dims(){[2]} * self.params[2].range,

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

    transpose := self >> TIterVStack(self.params[1].transpose(), self.params[2])
          .withTags(self.getTags()),

    isReal := self >> self.params[1].isReal(),

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

#    doNotMeasure := true,

    HashId := self >> let(
    	p := self.params,
    	h := When(IsBound(p[1].HashId), p[1].HashId(), p[1]),
            [h, p[2].range] :: When(IsBound(self.tags), self.tags, [])
    )
));

Class(TIterVStack, Tagged_tSPL_Container, rec(
    abbrevs :=  [ (nt, idx) -> Checked(IsSPL(nt), IsVar(idx),
	    [nt, idx]) ],

    dims := self >> (self.params[1].dims(){[1]} * self.params[2].range) :: self.params[1].dims(){[2]},

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

    transpose := self >> TIterHStack(self.params[1].transpose(), self.params[2])
          .withTags(self.getTags()),

    isReal := self >> self.params[1].isReal(),

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

#    doNotMeasure := true,

    HashId := self >> let(
    	p := self.params,
    	h := When(IsBound(p[1].HashId), p[1].HashId(), p[1]),
            [h, p[2].range] :: When(IsBound(self.tags), self.tags, [])
    )
));

Class(TNoDiagPullin, Tagged_tSPL_Container, rec(
    abbrevs :=  [ s -> Checked(IsSPL(s), [s]) ],
    HashId := self >> let(
    	p := self.params[1],
    	h := When(IsBound(p.HashId), p.HashId(), p),
            [ self.__name__, h ] :: When(IsBound(self.tags), self.tags, [])
        )
));

Class(TNoDiagPullinRight, Tagged_tSPL_Container, rec(
    abbrevs :=  [ s -> Checked(IsSPL(s), [s]) ],
    HashId := self >> let(
    	p := self.params[1],
    	h := When(IsBound(p.HashId), p.HashId(), p),
            [ self.__name__, h ] :: When(IsBound(self.tags), self.tags, [])
        )
));

Class(TNoDiagPullinLeft, Tagged_tSPL_Container, rec(
    abbrevs :=  [ s -> Checked(IsSPL(s), [s]) ],
    HashId := self >> let(
    	p := self.params[1],
    	h := When(IsBound(p.HashId), p.HashId(), p),
            [ self.__name__, h ] :: When(IsBound(self.tags), self.tags, [])
        )
));

Class(TNoPullRight, Tagged_tSPL_Container, rec(
    abbrevs :=  [ s -> Checked(IsSPL(s), [s]) ],
    HashId := self >> let(
    	p := self.params[1],
    	h := When(IsBound(p.HashId), p.HashId(), p),
            [ self.__name__, h ] :: When(IsBound(self.tags), self.tags, [])
        )
));

Class(TNoPullLeft, Tagged_tSPL_Container, rec(
    abbrevs :=  [ s -> Checked(IsSPL(s), [s]) ],
    HashId := self >> let(
    	p := self.params[1],
    	h := When(IsBound(p.HashId), p.HashId(), p),
            [ self.__name__, h ] :: When(IsBound(self.tags), self.tags, [])
        )
));

#F   tSPL ColMajor -> RowMajor transformation
Class(TColMajor, Tagged_tSPL_Container, rec(
    _short_print := true,
    abbrevs :=  [ (A) -> Checked(IsNonTerminal(A) or IsSPL(A), [A]) ],
    dims := self >> 2*self.params[1].dims(),
    terminate := self >> Mat(MatSPL(RC(self.params[1]))),

    transpose := self >> ObjId(self)(
	self.params[1].conjTranspose()).withTags(self.getTags()),

    conjTranspose := self >> self.transpose(),

    isReal := self >> true,

    # Do not use doNotMeasure, this will prevent TRC_By_Def from ever being found!
    doNotMeasure := false,
    normalizedArithCost := self >> self.params[1].normalizedArithCost(),

    HashId := self >> let(
	h := [ When(IsBound(self.params[1].HashId), self.params[1].HashId(),
		    self.params[1]) ],
        When(IsBound(self.tags), Concatenation(h, self.tags), h))
));


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

RulesFor(DCT4, rec(

    #F DCT4_Base2: DCT4_2 = (1,2) * R_13/8
    #F
    DCT4_Base2 := rec(
	info             := "DCT4_2 -> R_13/8",
	forTransposition := false,
	isApplicable     := P -> P[1] = 2,
	rule := (P, C) -> J(2) * Rot(13/8)
    ),

    #F DCT4_Base3: (base case for size 3)
    #F
    #F   DCT4_2 = (2,3) * (F_2 dirsum 1) * (a dirsum M) * (F_2 dirsum 1) * (2,3)
    #F
    #F derived by hand starting with AREP's decomposition by mon-mon symmetry.
    #F
    DCT4_Base3 := rec (
	info             := "DCT4_3 -> F_2, M, F_2",
	forTransposition := false,
	isApplicable     := P -> P[1] = 3,
	rule := (P, C) ->
	    Perm((2,3), 3) *
	    DirectSum(F(2), I(1)) *
	    DirectSum(Diag(Sqrt(3/8)), Sqrt(1/2)*Mat([[1/2, 1], [1, -1]])) *
	    DirectSum(F(2), I(1)) *
	    Perm((2,3), 3)
    ),

    #F DCT4_DCT2: DCT4_n = sums * DCT2_n * diag
    #F
    #F   inverse - transpose of RuleDCT4_2
    #F
    DCT4_DCT2 := rec (
	info             := "DCT4_n --> DCT2'_n",
	isApplicable     := P -> P[1] > 2,
	allChildren := P -> [[ DCT2(P[1]) ]],
	rule := (P, C) ->
	    sums1(P[1]) *
	    C[1] *
	    Diag(List([0..P[1]-1], i -> 1/(2*CosPi((2*i+1)/(4*P[1])))))
    ),

    #F DCT4_DCT2t: 1985, DCT4_n = sums * DCT2_n * diag
    #F
    #F   Chan:
    #F     Direct Methods for computing discrete sinusoidal transforms
    #F     IEE Proceedings, Vol. 137, 1990, pp. 433--442
    #F
    # switched off since critical path too long (sums6)
    # also there seems to be one mult too much
    DCT4_DCT2t := rec (
	info             := "DCT4_n --> DCT2'_n",
	switch           := false,
	isApplicable     := P -> P[1] > 2,
	allChildren      := P -> [[ DCT2(P[1]) ]],
	rule := ( P, C ) ->
	    sums6(P[1]) *
	    C[1] *
	    Diag(List([0..P[1] - 1], i -> 2 * CosPi((2*i + 1)/(4*P[1]))))
    ),

    #F DCT4_DCT2andDST2: 1985
    #F
    #F   DCT4_n = blocks * perm * (DCT2_n/2 dirsum DST2_n/2) * rotations
    #F
    #F   Wang: (transposed)
    #F     On Computing the Discrete Fourier ans Cosine Transforms,
    #F     IEEE Trans. on Signal Proc., 1985, pp. 1341--1344
    #F
    DCT4_DCT2andDST2 := rec (
	info             := "DCT4_n --> DCT2'_n/2, DST2'_n/2",
	isApplicable     := P -> IsInt(P[1]) and P[1] > 2 and P[1] mod 2 = 0,
	allChildren := P -> [[ DCT2(P[1]/2), DST2(P[1]/2) ]],
	rule := (P, C) ->
	    DirectSum(I(1),
		      Tensor(I(P[1]/2-1), F(2)*J(2)),
		      I(1)) *
	    DirectSum(C[1], C[2]) ^ L(P[1], 2) *
	    DirectSum(List([0..P[1]/2-1], i -> J(2)*Rot(2-(2*P[1]-2*i-1)/(4*P[1])))) *
	    LIJ(P[1])

	  # LIJ(P[1]) *
	  # DirectSum(C[1], C[2] ^ J(P[1]/2)) *
	  # DirectSum(List([0..P[1]/2-1], i -> J(2)*Rot(2-(2*P[1]-2*i-1)/(4*P[1]))))
	  #    ^ LIJ(P[1])
    ),

    #F DCT4_DST4andDST2: 1988
    #F
    #F   DCT4_n = diag * sums * perm * (DST4_n/2 dirsum DST2_n/2) * sign *
    #F            (DFT_2 tensor I_n/2) * diag * perm
    #F
    #F   Rao/Yip:
    #F     The Decimation-in-Frequency Algorithms for a family of
    #F     Discrete Sine and Cosine Algorithms.
    #F     Circuits, Systems, and Signal Processing, 1988, pp. 3--19
    #F
    DCT4_DST4andDST2 := rec (
	info             := "DCT4_n --> DST4_n/2, DST2_n/2",
	switch           := false,
	isApplicable     := P -> IsInt(P[1]) and P[1] > 2 and P[1] mod 2 = 0,

	allChildren := P -> [[ DST4(P[1]/2), DST2(P[1]/2) ]],

	rule := (P, C) ->
	    Diag(List([0..P[1] - 1], i -> (-1)^i)) *
	    sums1(P[1]).transpose() *
	    L(P[1], P[1]/2) *
	    DirectSum(C[1], C[2] * (-1)) *
	    Tensor(F(2), I(P[1]/2)) *
	    Diag(Concat(List([1..P[1]/2], i -> 1/(2*SinPi((2*i - 1)/(4*P[1])))),
		        List([1..P[1]/2], i -> 1/(2*CosPi((2*i - 1)/(4*P[1])))))) *
	    IJ(P[1], P[1]/2)
    ),

    #F DCT4_Iterative: 1985, DCT4_n = iterative
    #F
    #F   Chen/Smith/Fralick:
    #F     A Fast Computational Algorithm for the Discrete Cosine Transform,
    #F     IEEE Trans. on Comm., 1977, pp. 1004--1009
    #F   corrected in:
    #F   Wang:
    #F     Reconsideration of --above--, IEEE Trans. on Comm., 1983, pp. 121--123
    #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
    DCT4_Iterative := rec (
	info             := "DCT4_n iterative",
	isApplicable     := P -> P[1] > 2 and Is2Power(P[1]),
	rule := (P, C) ->
	    LIJ(P[1]) * L(P[1], 2) *
	    DirectSum(List([1..P[1]/2], i -> Rot(1/2-(4*i-3)/(4*P[1])) * J(2))) *

	    Compose(List([1..Log2Int(P[1])-1], j -> let(jj := 2^j,
		    Tensor(I(jj/2), F(2), I(P[1]/jj)) *
		    Tensor(I(jj/2),
			DirectSum(I(P[1]/jj), List([1..P[1]/jj/2],
			     i -> Rot(1/2-(4*i-3)/(2*P[1]/jj)) * J(2))))))) *
	    perm6(P[1])
    )
));


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

SetLinuxX86Profile:=function(opts, p, use_openmp)
   opts.profile:=Copy(default_profiles.linux_x86_icc);
   SetMakeOptsAffinity(opts);
   SetMakeOptsLibgen(opts);
   SetMakeOptsAssembly(opts);
   opts.includes := Concat(opts.includes, [ "<string.h>"]);

   if IsBound(opts.vector) then
       SetMakeOptsSSE(opts);
   fi;

   if p>1 then
       if (use_openmp) then
           SetMakeOptsOpenMP(opts);
       else
           SetMakeOptsPThreads(opts);
       fi;
   fi;
end;

SetOLOptions:=function(isa, p, use_openmp)
   local opts, MixedUnparser;

   if (isa = false) or (IsString(isa)) then
       opts := Copy(SpiralDefaults);
       if (IsString(isa)) then opts := InitDataType(opts,isa); fi;
   else
       if (isa = true) then isa:=SIMD_ISA_DB.active()[1]; fi;
       opts := SIMDGlobals.getOpts(isa);
   fi;

   if not IsBound(opts.profile) then
       SetLinuxX86Profile(opts, p, use_openmp);
   else
#Hack Cell profile
       opts.profile := default_profiles.linux_cellSPU_gcc_MMM;
       opts.profile.stubopts.ROWS := 1;
       opts.profile.stubopts.COLUMNS := 1;
       opts.includes := Concat(opts.includes, [ "<mm_malloc.h>", "<string.h>"]);
       opts.includes := Filtered(opts.includes, x->x<>"<omega32.h>");
   fi;

   if p>1 then
       if (use_openmp) then
           opts.unparser:=Class(MixedUnparser,OpenMP_UnparseMixin,opts.unparser);
       else
           opts.unparser:=Class(MixedUnparser,SMP_UnparseMixin,opts.unparser);
           opts.subParams := [var("num_threads", TInt), var("tid", TInt)];
       fi;
   fi;

   #OL specific
   opts.codegen := OLCodegen;

   opts.formulaStrategies.sigmaSpl := OLDefaultStrategy;

   opts.formulaStrategies.postProcess := [
              OLVectorPropagateRuleset,
              OLPushScatQuestionMarkInRuleset,
              OLAlreadyInitializedScatQuestionMarkRuleset,
              OLCrossPullInRuleset,
              OLAfterCrossPullInRuleset,
              OLVectorPropagateRuleset,
              OLSingleComposeRuleset,
              (s, opts) -> compiler.BlockSumsOpts(s, opts)
          ];

   #Unrolling
   opts.markBlock := MarkBlocksOps;
   opts.globalUnrolling := 300;

   #compile options
   opts.compileStrategy :=IndicesCS2;
   opts.useDeref := true;
   opts.propagateNth := false;
   opts.doScalarReplacement := true;

   #Libgen requires full spec of perms
   VPerm.print:=VPerm.printl;

   #ScatAcc requires full zero allocation for now
   opts.zeroallocate := true;

   #final code options
   opts.subName:="multi";
   opts.subInitName:="init_multi";

   return opts;
end;

SetSAROptions:=function(isa, p, use_openmp)
  local tags,opts;
  opts := SetOLOptions(isa, p, use_openmp);
  opts.codegen := OLCodegen;
  opts.zeroallocate := false;
  opts.InputTypes := [TDouble,TDouble,TDouble,TDouble];
  opts.OutputTypes := [TDouble,TDouble];
  opts.generateComplexCode := true;
  opts.TRealCType := "complex";
  opts.TRealCtype := "complex";
  opts := InitDataType(opts,"f64c");
  opts.formulaStrategies.sigmaSpl := OLDefaultStrategy;
  opts.breakdownRules.TTensorI_OL := [TTensorI_OL_Base, TTensorI_OL_Parrallelize_AParFirst, TTensorI_OL_Vectorize_AVecLast];
  opts.breakdownRules.SAR := [SAR_Base,SAR_MatchFilter,SAR_MatchFilterSMP,SAR_MatchFilterVec,SAR_MatchFilterVecBase];
  opts.breakdownRules.ExpandMult := [ExpandMult_Base, ExpandMult_One];
  opts.breakdownRules.SAR_Interpolation := [SAR_Interpolation_One];
  tags:=[];

  if (p>1) then
      Add(tags,AParSMP(p));
  fi;
  opts.formulaStrategies.postProcess := [
        OLVectorPropagateRuleset,
 	      OLScatQuestionMarkToScat,
              OLPushScatAccRuleset,
              OLCrossPullInRuleset,
              OLAfterCrossPullInRuleset,

              OLVectorPropagateRuleset,

              OLSingleComposeRuleset,
	      (s, opts) -> compiler.BlockSumsOpts(s, opts)];

  #Unrolling
  opts.markBlock := MarkBlocksOps;
  opts.globalUnrolling := 300;

  #compile options
  opts.compileStrategy :=IndicesCS2;
  opts.useDeref := true;
  opts.doNotScalarize := false;
  opts.propagateNth := false;
  opts.doScalarReplacement := true;

  #ScatAcc requires full zero allocation for now
  opts.zeroallocate := false;

  #final code options
  opts.subName:="multi";
  opts.subInitName:="init_multi";

  return opts;

end;


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

Declare(VStack);

# ==========================================================================
# _HStack(<spl1>, ...)
#    equivalent of HStack, expressed using Tensor/Mat/DirectSum
#    use for unrolled code, to avoid problems with array scalarization
#
_HStack := arg -> Tensor(Mat([Replicate(Length(arg), 1)]), I(Rows(arg[1]))) * DirectSum(arg);

# ==========================================================================
# _SUM(<spl1>, ...)
#    overlapping dense spl1+spl2+..., expressed using Tensor/Mat/VStack
#    use for unrolled code, to avoid problems with array scalarization
#
_SUM    := arg -> Tensor(Mat([Replicate(Length(arg), 1)]), I(Rows(arg[1]))) * VStack(arg);

# ==========================================================================
# HStack(<spls>) - horizontal stacking operator (row of blocks)
#
# Example: MatSPL( HStack(F(2), I(2))) );
#     [ [ 1, 1, 1, 0 ],
#       [ 1, -1, 0, 1 ] ]
# ==========================================================================
Class(HStack, BaseOperation, rec(
    abbrevs := [ arg -> [Flat(arg)] ],
    new := (self, spls) >> Checked(
		IsList(spls), Length(spls) >= 1,
		SPL(WithBases(self, rec(_children := spls)).setDims())
    ),

    dims := self >> [Maximum(List(self.children(), c->c.dims()[1])), Sum(self._children, x->x.dims()[2])],

    isPermutation := self >> false,
    transpose := self >> CopyFields(self, VStack(List(self._children, x->x.transpose()))),
    conjTranspose := self >> CopyFields(self, VStack(List(self._children, x->x.conjTranspose()))),

    toAMat := self >> ApplyFunc(_HStack, self.children()).toAMat()
));

# this is a HStack that is guaranteed not to perform any ops, and doesnt need ScatAcc in SumsGen
Class(HStack1, HStack);

# ==========================================================================
# VStack(<spls>) - vertical stacking operator (column of blocks)
#
# Example: MatSPL( VStack(F(2), I(2))) );
#    [ [ 1, 1 ],
#      [ 1, -1 ],
#      [ 1, 0 ],
#      [ 0, 1 ] ]
# ==========================================================================
Class(VStack, BaseOperation, rec(
    abbrevs := [ arg -> [Flat(arg)] ],

    new := (self, spls) >> Checked(
		IsList(spls), Length(spls) >= 1,
		SPL(WithBases(self, rec(_children := spls)).setDims())
    ),

    dims := self >> [Sum(self._children, x->x.dims()[1]), Maximum(List(self.children(), c->c.dims()[2]))],

    isPermutation := self >> false,
    transpose := self >> CopyFields(self, HStack(List(self._children, x->x.transpose()))),
    conjTranspose := self >> CopyFields(self, HStack(List(self._children, x->x.conjTranspose()))),

    toAMat := self >> ApplyFunc(_HStack, List(self.children(), e -> e.transpose())).transpose().toAMat()
));

# ==========================================================================
# BlockMat(<matrix of spls>) - puts several spls in a bigger block matrix
#
# Example: MatSPL( BlockMat([[I(2), J(2)], [O(2,2), -1 * I(2)]]) )
#   [ [ 1,  0,  0,  1 ],
#     [ 0,  1,  1,  0 ],
#     [ 0,  0, -1,  0 ],
#     [ 0,  0,  0, -1 ] ]
# ==========================================================================
BlockMat := arg -> let(mat := When(Length(arg)=1, arg[1], arg),
    VStack(List(mat, x->HStack(x))));