text:
########################################################################
#
# walrus: Computational Methods for Finitely Generated Monoids and Groups
#
# WARNING
#
########################################################################
#
# TODO (documentation)
#  - Examples
#  - notes on the implementation (design decisions etc)
#
# TODO (technical)
#  - MakeIsReducedDiagram an operation
#  - Sort out families of elements etc
#  - what to do about the __ID hack?
#  - What to do about float vs rational?
#  - provide a way to label generators
#  - Introduce an InfoLevel and log things with info levels, remove Print
#    statements
#  - Put relators of different presentations into different families
#  - Check whether some sub-functions could be cleaned out from functions
#    and would become more generally useful
#  - Better indexing of locations on relator/places
#  - Better datastructure for storing one step reachables
#    (basically a hashmap from orb I think) plus wrapper functions
#  - PlaceQuadruples should be an object or a record
#  - Add proper error messages for "This shouldn't happen" with a pointer
#    towards the github issue tracker
#  - Make a reasonably standard format (and exporter/importer) to exchange
#    pregroup presentations so we can move them between GAP/other
#    implementations for cross-checking

#
# TODO (mathematical/functional)
#  - Write tests
#  - Interleaving
#  - Preparing the presentation (User choosing Pregroup/removing generators, relators)
#  - check LocationBlobGraph + distances
#  - Make the distance computation in LBG more efficient?
#  - Check OneStepReachables
#  - Check RSymTest
#
########################################################################
#
# Some notes on the implementation
#
# - Stored curvature is always "negative curvature", so in most cases a
#   positive value.
# - Curvature values are stored as rationals at first, later converted to
#   floats. This might not be the best way to do this (we might want to
#   always use rationals or always use floats)
#

# Alternative to LocationBlobGraph, only has pairs, coloured Green for location pair
# and red for intermult pair
# a bit hacky though
InstallMethod(VertexGraph, "for a pregroup presentation",
              [IsPregroupPresentation and IsPregroupPresentationRep],
function(pres)
    local mat, v, vd, vg, e, loc, loc2, pg, p, ploc, ploc2, imp, i, backmap, len, li, labels;

    pg := Pregroup(pres);

    li := LocationIndex(pres);

    # Green vertices [[a, b], 0] correspond to locations
    v := List(Keys(li), l -> [l, 0]);
    # Red vertices correspond to intermult pairs
    Append(v, List(IntermultPairs(pg), p -> [p, 1]));

    e := function(a,b)
        if (a[2] = 0) and (b[2] = 0) then              # both vertices are green
            if PregroupInverse(a[1][2]) = b[1][1] then # There are locations R(i,a,b) and R'(j,b^-1,c)
                ploc := li[a[1]];
                ploc2 := li[b[1]];
                for loc in ploc do
                    for loc2 in ploc2 do
                        if CheckReducedDiagram(loc, loc2) then
                            return 1;
                        fi;
                    od;
                od;
                return infinity;
            fi;
        elif (a[2] = 0) and (b[2] = 1) then            # a is green, b is red
            if PregroupInverse(a[1][2]) = b[1][1] then # There is location R(i,a,b)
                return 1;
            fi;
        elif (a[2] = 1) and (b[2] = 0) then            # a is red, b is green
            if PregroupInverse(a[1][2]) = b[1][1] then #
                return 0;
            fi;
        fi;
        return infinity;
    end;

    backmap := HashMap();
    for i in [1..Length(v)] do
        backmap[ [__ID(v[i][1][1]), __ID(v[i][1][2]), v[i][2] ] ] := i;
    od;

    vg := Digraph(v, {i,j} -> e(i,j) <> infinity);
    SetDigraphVertexLabels(vg, v);

    mat := List(DigraphVertices(vg), i -> List(DigraphVertices(vg), j -> e(v[i], v[j]) ) );
    vg!.mat := mat;
    vg!.backmap := backmap;

    return vg;
end);

InstallMethod(VertexGraphDistances, "for a pregroup presentation",
              [IsPregroupPresentation and IsPregroupPresentationRep],
function(pres)
    local vg, wt, f, vertices, n, mat, i, j, k, a, b, t, out;

    vg := VertexGraph(pres);
    vertices := DigraphVertices(vg);
    mat := vg!.mat;

#    SC_FLOYD_WARSHALL(mat);

    for k in vertices do
        for i in vertices do
            for j in vertices do
                t := mat[i, k] + mat[k, j];
                if t < mat[i,j] then
                    mat[i,j] := t;
                fi;
            od;
        od;
    od;
    # Can there be more than one such entry?
    return mat;
end);

InstallMethod(VertexTriples, "for a pregroup presentation",
            [IsPregroupPresentation and IsPregroupPresentationRep],
function(pres)
    local v, v1, v2, v1l, v2l, vl, vg, lv, vgd, dist, vtl;

    vg := VertexGraph(pres);
    vgd := VertexGraphDistances(pres);
    vtl := HashMap(16384);

    for v in DigraphVertices(vg) do

        vl := DigraphVertexLabel(vg, v);
        # Only green vertices
        if (vl[2] = 0) then
          for v1 in InNeighboursOfVertex(vg, v) do
              v1l := DigraphVertexLabel(vg, v1);
              for v2 in OutNeighboursOfVertex(vg, v) do
                  v2l := DigraphVertexLabel(vg, v2);
                  dist := vgd[v2][v1];
                  if (v1l[2] = 0) and (v2l[2] = 0) then
                      if dist = 1 then
                          vtl[ [v, v1, v2 ] ] := 1/6;
                      elif dist = 2 then
                          vtl[ [ v, v1, v2 ] ] := 1/4;
                      elif dist = 3 then
                          vtl[ [ v, v1, v2 ] ] := 3/10;
                      elif dist > 3 then
                          vtl[ [ v, v1, v2 ] ] := 1/3;
                      else
                          Error("this shouldn't happen");
                      fi;
                  elif (v1l[2] = 0) and (v2l[2] = 1) then
                      if dist = 0 then
                          vtl[ [ v, v1, v2 ] ] := 0;
                      elif dist = 1 then
                          vtl[ [ v, v1, v2 ] ] := 1/6;
                      elif dist > 1 then
                          vtl[ [ v, v1, v2 ] ] := 1/4;
                      else
                          Error("this shouldn't happen");
                      fi;
                  elif (v1l[2] = 1) and (v2l[2] = 0) then
                      if dist = 1 then
                          vtl[ [ v, v1, v2 ] ] := 0;
                      elif dist = 2 then
                          vtl[ [ v, v1, v2 ] ] := 1/6;
                      elif dist > 2 then
                          vtl[ [ v, v1, v2 ] ] := 1/4;
                      else
                          Error("this shouldn't happen");
                      fi;
                  elif (v1l[2] = 1) and (v2l[2] = 1) then
                      vtl[ [ v, v1, v2 ] ] := 0;
                      # what if dist=infinity (i.e. no path)
                  else
                      Error("this should not happen");
                  fi;
              od;
          od;
        fi;
    od;
    return vtl;
end);

InstallGlobalFunction(Vertex,
function(pres, v1, v, v2)
    local vt, t;

    t := VertexTriples(pres)[ [v, v1, v2] ];
    if t <> fail then
        return t;
    else
        return 1/3;
    fi;
end);

#XXX Computes triples (a,b,c) that are infixes
#    of strings found here with appropriate numbers
InstallMethod(ShortRedBlobIndex, "for a pregroup presentation",
              [IsPregroupPresentation],
function(pres)
    local pg, n, imm, index, nrl, cand, word, nonrletts, len, c;

    pg := Pregroup(pres);
    n := Size(pg);
    imm := IntermultMapIDs(pg);

    # TODO: This would be better done by a tree structure
    index := HashMap();

    cand := [[2..n]];

    nonrletts := [];
    word := [];

    len := 1;

    # minimum length 3
    while (len > 0) and (len < 7) do
        if Length(cand[len]) > 0 then
            c := Remove(cand[len], 1);
            word[len] := pg[c];
            nonrletts[len] := not IsRLetter(pres, pg[c]);
            # Possible not go further if we have too many R-letters?
            if IsIntermultPair(word[len], word[1]) and
               ReduceUPregroupWord(word{ [1..len] }) = [] then
                nrl := SizeBlist(nonrletts{[1..len]});
                if len = 3 then
                    # if len = 3 then no proper subword can be equal to 1
                    # because we only look at intermult pairs, and pairs
                    # of inverses are explicitly not intermult pairs
                    if nrl = 0 then
                        EnterAllSubwords(index, word{[1..len]}, 1/6);
                    elif nrl = 1 then
                        EnterAllSubwords(index, word{[1..len]}, 1/4);
                    else
                        # Print("? nrl", nrl, "\n");
                    fi;
                elif len = 4 then
                    if nrl = 0 then
                        EnterAllSubwords(index, word{[1..len]}, 1/4);
                    elif nrl = 1 then
                        EnterAllSubwords(index, word{[1..len]}, 1/3);
                    else
                        # Print("? nrl", nrl, "\n");
                    fi;
                elif len = 5 then
                    if (nrl = 0) then
                        EnterAllSubwords(index, word{[1..len]}, 3/10);
                    fi;
                elif (len = 6)
                     and (ReduceUPregroupWord(word{[2,3,4]}) <> [])
                     and (ReduceUPregroupWord(word{[3,4,5]}) <> []) then
                    if nrl = 0 then
                        EnterAllSubwords(index, word{[1..len]}, 1/3);
                    else
                        # Print("? nrl", nrl, "\n");
                    fi;
                fi;
                # We have entered, so backtrack
                len := len - 1;
            else
                if (len = 6) then
                    len := len - 1;
                else
                    len := len + 1;
                    cand[len] := ShallowCopy(imm[c]);
                fi;
            fi;
        else
            len := len - 1;
        fi;
    od;
    return index;
end);

# This should become an operation
InstallMethod(Blob, "for a pregroup presentation, an element, an element, and an element"
              , [ IsPregroupPresentation
                , IsElementOfPregroup
                , IsElementOfPregroup
                , IsElementOfPregroup ],
function(pres, a, b, c)
    local it, v, n;

    n := Size(Pregroup(pres));
    it := ShortRedBlobIndex(pres);
    v := it[ [__ID(a), __ID(b), __ID(c)] ];

    if v <> fail then
        return v;
    fi;
    return 5/14;
end);

# TODO: Use a map
InstallGlobalFunction(VertexFor,
function(vg, trip)
    return vg!.backmap[ [ __ID(trip[1]), __ID(trip[2]), trip[3] ] ];
end);

# Compute a list of
#  - places reachable from loc1 on Relator(loc1), and
#  - the corresponding location on Relator(loc2)
# by a consolidated edge between Relator(loc1) and Relator(loc2)
InstallGlobalFunction(ConsolidatedEdgePlaces,
function(loc1, loc2)
    local res, length, r1_loc, r2_loc, r1_length;

    res := [];
    length := 0;

    r1_loc := loc1;
    r2_loc := loc2;
    r1_length := Length(Relator(loc1));

    if OutLetter(r1_loc) <> PregroupInverse(InLetter(r2_loc)) then
        Error("This shouldn't happen");
#        return [];
    fi;

    repeat
        r1_loc := NextLocation(r1_loc);
        r2_loc := PrevLocation(r2_loc);
        length := length + 1;
        Append(res, List(Places(r1_loc), x -> [x, r2_loc, length]));
    until (OutLetter(r1_loc) <> PregroupInverse(InLetter(r2_loc)))
          or (length = r1_length);

    # Meh.
    if (length = r1_length) or
       (length = Length(Relator(loc2))) then
        return [];
    else
        return res;
    fi;
end);

# Helper function to add or update
# a "place quadruple"
AddOrUpdatePQ := function(L, pq)
    local pp;

    pp := PositionProperty(L, x -> x[1] = pq[1] and x[2] = pq[2]);
    if pp = fail then
        Add(L, Immutable(pq));
    else
        if L[pp][4] < pq[4] then
            L[pp] := Immutable(pq);
        fi;
    fi;
end;

# The RSym tester
# The epsilon is chosen by the user
InstallMethod(RSymTestOp, "for a pregroup presentation, and a rational",
              [IsPregroupPresentation, IsRat],
function(pres, eps)
    local i, j, rel, R, pplaces,
          places, Ps, P, Q, Pq,
          osrp, # a one-step reachable place
          osrpp,
          L,
          zeta,
          xi, psip, pp;
    # Make sure epsilon is a float
    # FIXME: To experiment and find out about rounding errors
    #        we try running this thing with rationals
    eps := Rat(eps);
    # FIXME: CAREFUL, if the zeta value is wrong, the RSymTester
    #        does not work correctly
    zeta := Minimum(Int(6 * (1 + eps) + 1) - 1,
                    LengthLongestRelator(pres));
    Info(InfoWalrus, 10
         , "RSymTest start");
    Info(InfoWalrus, 10
         , "zeta: ", zeta);

    # Precompute some things so it doesn't show up in weird places
    # in the profiling data
    # Some of this we might want to compute on demand?
    VertexGraphDistances(pres);
    pplaces := Places(pres);
    if not IsEmpty(pplaces) then
        OneStepReachablePlaces(pplaces[1]);
    fi;

    for rel in Relators(pres) do
        Info(InfoWalrus, 20
             , "relator: ", ViewString(rel), "\n");
        places := Places(rel);
        for Ps in places do
            Info(InfoWalrus, 20
                 , "  start place: "
                 , ViewString(Ps)
                 , "\n" );
            L := [ [__ID(Ps), 0, 0, 0] ]; # This list is called L in the paper
            # This is the list of possible decompositions
            # The meaning of the components of the quadruples q is
            # - q[1] is a place
            # - q[2] is the distance of q[1] from Ps
            # - q[3] is the number of steps that q[1] is from Ps
            # - q[4] is a curvature value
            for i in [1..zeta] do
                Info(InfoWalrus, 30
                     , STRINGIFY("L = ", ViewString(L)), ", i = ", i);
                for Pq in L do      # Pq is for "PlaceQuadruple", which is
                                    # a silly name
                    Info(InfoWalrus, 30,
                         STRINGIFY("Considering: ", Pq));
                    if Pq[3] = i - 1 then  # Reachable in i - 1 steps
                        osrpp := OneStepReachablePlaces(pplaces[Pq[1]]);
                        for osrp in Keys(osrpp) do
                            Info(InfoWalrus, 30,
                                 STRINGIFY("OneStepReachable: ", osrp));
                            if Pq[2] + osrp[2] <= Length(rel) then
                                psip := Pq[4]
                                        + osrp[2] * (1 + eps) / Length(rel)
                                        # storing positive values -> subtract here
                                        - osrpp[osrp];
                                Info(InfoWalrus, 30
                                     , STRINGIFY("psi' = "
                                                , Pq[4], " + "
                                                , osrp[2] * (1+eps) / Length(rel), " - "
                                                , osrpp[osrp], " = "
                                                , psip, "\n")
                                    );
                                if psip < 0 then
                                elif (Pq[4] > 0) and
                                     (osrp[1] = __ID(Ps)) and
                                     (Pq[2] + osrp[2] = Length(rel)) then
                                    # Add the place that caused the failure
                                    Add( L, Immutable( [ osrp[1], Pq[2] + osrp[2], i, psip ] ) );
                                    return [fail, L, Pq];
                                else
                                    AddOrUpdatePQ(L, [osrp[1], Pq[2] + osrp[2], i, psip]);
                                fi;
                            fi;
                        od;
                    fi;
                od;
            od;
        od;
    od;
    return true;
end);

InstallGlobalFunction(RSymTest,
function(args...)
    local pg, pgp, rgreen, inv;

    if Length(args) >= 2 then
        if IsPregroupPresentation(args[1]) then
            return RSymTestOp(args[1], args[2]);
        fi;
        if IsFreeGroup(args[1]) then
            pgp := PregroupPresentationFromFp(args[1], args[2], args[3]);
            return RSymTestOp(pgp, args[4]);
        fi;
    else
        # TODO: Usage message
    fi;
end);

InstallMethod(IsHyperbolic, "for a free group, a list, a list, and a rational number",
              [IsFreeGroup, IsObject, IsObject, IsRat],
              {freeg, rrel, grel, eps} -> RSymTest(freeg, rrel, grel, eps));

InstallMethod(IsHyperbolic, "for a pregroup presentation, and a rational number",
              [IsPregroupPresentation, IsRat],
              {pres, eps} -> RSymTest(pres,eps));

ext: gi

max_stars_repo_name: wilfwilson/walrus

max_issues_repo_name: wilfwilson/walrus

max_forks_repo_name: wilfwilson/walrus


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

Class(Drawer, rec(
    __call__ := meth(self,n)
        return WithBases(self,rec(n:=n,instances:=[]));
    end,

    add := meth(self,instance)
        if (not IsList(instance)) or (Length(instance)<>self.n) then
            Error("Trying to add a non regular instance");
        fi;
        Add(self.instances,instance);
    end,

    getNormalTransfo := meth(self)
       local result,l,m,s;
       result:=[];
       for l in Transposed(self.instances) do
          m:=Mean(l);
          s:=Try(StdDev(l));
          if not s[1] then s:=0; else s:=s[2]; fi;
          Add(result,[m,s]);
       od;
       return result;
    end,

    normalize := meth(self,transfo)
       local t,n,l;
       t:=[];
       n:=1;
       for l in Transposed(self.instances) do
          if (transfo[n][2]<>0) then
              Add(t,List(l,x->(x-transfo[n][1])/transfo[n][2]));
          else
              Add(t,List(l,x->(x-transfo[n][1])));
          fi;
          n:=n+1;
       od;
       self.instances:=Transposed(t);
    end,

    checkCredibility := meth(self,testing, transfo)
        local j;
        for j in [1..Length(self.instances)] do
            for i in [1..Length(testing.instances[1])] do
               if transfo[i][2]=0 and transfo[i][1]<>testing.instances[j][i] then
                   self.instances[j]:=[2^16,1,2^32];
               fi;
            od;
        od;
    end,

    denormalize := meth(self,transfo)
       local t,n,l;
       t:=[];
       n:=1;
       for l in Transposed(self.instances) do
          if (transfo[n][2]<>0) then
              Add(t,List(l,x->x*transfo[n][2] + transfo[n][1]));
          else
              Add(t,List(l,x->x + transfo[n][1]));
          fi;
          n:=n+1;
       od;
       self.instances:=Transposed(t);
    end,

    print := meth(self,filename)
       local l,a;
       PrintTo(filename,"");
       for l in self.instances do
          for a in l do
             AppendTo(filename,a," ");
          od;
          AppendTo(filename,"\n");
       od;
    end,

));

Class(GaussianPredictor, rec(

    __call__ := meth(self)
        return WithBases(self);
    end,

    predict := meth(self, hash, objid, testCases, tags)
      local n,a,i,transfo,training,testing,predict,index,sol,t;
      n:=DimensionsMat(testCases)[2];
      training:=Drawer(n+1);
      for a in Flat(hash.entries) do
          if ObjId(a.key)=objid and a.key.params[n+1]=tags and Length(a.data)>0 and IsBound(a.data[1].mflopslibgen) then
              training.add(Concat(a.key.params{[1..n]},
                      [a.data[1].mflopslibgen]));
          fi;
      od;

      transfo:=training.getNormalTransfo();
      a:=Copy(training);
      a.normalize(transfo);
      a.print("trainingset");

      testing:=Drawer(n);
      testing.instances:=testCases;
      a:=Copy(testing);
      a.normalize(DropLast(transfo,1));
      a.print("testingset");

      IntExec("rm -f rlog.gp");
      IntExec(Concat("gp-spec rlog.gp ",String(n),
              " 1 10 / 0.05:0.4 0.05:0.4"));
      IntExec("model-spec rlog.gp real 0.05:0.4");
      IntExec(Concat("data-spec rlog.gp ",String(n)," 1 / trainingset ."));
      IntExec("gp-gen rlog.gp fix 0.5 0.1");
      IntExec("mc-spec rlog.gp heatbath hybrid 18:3 0.5");
      IntExec("gp-mc rlog.gp 100");
      IntExec("gp-pred nqb rlog.gp 51: / testingset >myoutput");
      IntExec("awk 'BEGIN {print \"[\"} {print \"[\",$1,\", \",$2,\", \",$3,\"], \"} END{print\"];\"}' myoutput >myrealoutput");

      a:=ReadVal("myrealoutput");

      predict:=Drawer(3);
      predict.instances:=a;
      predict.denormalize(Replicate(3,Last(transfo)));

      for a in [1..Length(testing.instances)] do
          t:=ApplyFunc(objid,Concat(testing.instances[a],[tags]));
          i:=HashLookup(hash,t);
          if (i<>false) and IsBound(i[1].mflopslibgen) then
              predict.instances[a]:=Replicate(3,i[1].mflopslibgen);
          fi;
      od;

      return predict.instances;
   end,

    predict2 := meth(self, hash, objid, testCases, tags)
      local n,a,i,transfo,training,testing,predict,index,sol,t,tmp;
      n:=DimensionsMat(testCases)[2];
      training:=Drawer(n+1);
      for a in Flat(hash.entries) do
          if ObjId(a.key)=objid and a.key.params[n+1]=tags and Length(a.data)>0 and
              IsBound(a.data[1].measured)  then
              training.add(Concat(a.key.params{[1..n]},
                      [_compute_mflops(When(IsBound(a.key.normalizedArithCost),
                                  a.key.normalizedArithCost(), 0), a.data[1].measured)]));
          fi;
      od;
      Print("Estimating ",Length(testCases)," points (with ",Length(training.instances),"): ",
          objid," (",tags,")\n");

      if Length(training.instances)=0 then
          return Replicate(Length(testCases),[2^16,0,2^32]);
      fi;
      transfo:=training.getNormalTransfo();
      a:=Copy(training);
      a.normalize(transfo);
      a.print("trainingset");

      testing:=Drawer(n);
      testing.instances:=testCases;
      a:=Copy(testing);
      a.normalize(DropLast(transfo,1));
      a.print("testingset");

      IntExec("rm -f rlog.gp");
      IntExec(Concat("gp-spec rlog.gp ",String(n),
              " 1 10 / 0.05:0.4 0.05:0.4"));
      IntExec("model-spec rlog.gp real 0.05:0.4");
      IntExec(Concat("data-spec rlog.gp ",String(n)," 1 / trainingset . 2>/dev/null"));
      IntExec("gp-gen rlog.gp fix 0.5 0.1");
      IntExec("mc-spec rlog.gp heatbath hybrid 18:3 0.5");
      IntExec("gp-mc rlog.gp 100");
      IntExec("gp-pred nqb rlog.gp 51: / testingset >myoutput");
      IntExec("awk 'BEGIN {print \"[\"} {print \"[\",$1,\", \",$2,\", \",$3,\"], \"} END{print\"];\"}' myoutput >myrealoutput");

      a:=ReadVal("myrealoutput");

      predict:=Drawer(3);
      predict.instances:=a;
      predict.denormalize(Replicate(3,Last(transfo)));

      predict.checkCredibility(testing, transfo);

      for a in [1..Length(testing.instances)] do
          t:=ApplyFunc(objid,Concat(testing.instances[a],[tags]));
          i:=HashLookup(hash,t);
          if (i<>false) and IsBound(i[1].measured) then
              predict.instances[a]:=Replicate(3,i[1].measured);
          else
              predict.instances[a]:=List(predict.instances[a],x->
                  _compute_cycles_from_mflops(When(IsBound(t.normalizedArithCost),
                                  t.normalizedArithCost(), 0),x));
              #invert max and min because min flops is max runtime
              tmp:=predict.instances[a][2];
              predict.instances[a][2]:=predict.instances[a][3];
              predict.instances[a][3]:=When(tmp<0,2^64,tmp);
          fi;
      od;

      return predict.instances;
   end,

   countTraining:=meth(self, dpbench, objid, n, tags)
      local r,hash,a;
      hash:=dpbench.exp.bases.hashTable;
      r:=0;
      for a in Flat(hash.entries) do
          if ObjId(a.key)=objid and a.key.params[n+1]=tags and Length(a.data)>0 then
              r:=r+1;
          fi;
       od;
      return r;
   end,

   findMax:=meth(self, base_bench, libgen_bench, base_objid, libgen_objid, testing, tags,flopsfactor)
       local over,sol,secured,estimate,a,s, unscaledsol;
       base_bench.verbosity:=0;
       if (g.countTraining(base_bench,base_objid,DimensionsMat(testing)[2],tags)<2) then
           base_bench.exp.bases.benchTransforms:=[ApplyFunc(base_objid,Concat(First(testing,x->true),[tags])),[ApplyFunc(base_objid,Concat(Last(testing),[tags]))]];
           base_bench.resumeAll();
       fi;

       over:=false;
       while(not over) do
           unscaledsol:=self.predict(base_bench.exp.bases.hashTable,base_objid,testing,tags);
           sol:=List([1..Length(testing)],x->List(unscaledsol[x],y->flopsfactor(testing[x],y)));
           secured:=[testing[1],sol[1]];
           estimate:=[testing[1],sol[1]];
           for a in [1..Length(testing)] do
              s:=[testing[a],sol[a]];
              if (s[2][1]=s[2][2] and s[2][1]=s[2][3] and s[2][1]>secured[2][1]) then
                  secured:=s;
              fi;
              if (s[2][3]>estimate[2][3]) then
                  estimate:=s;
              fi;
           od;

           Print("Current Maximum is ",Red(StringDouble("%.3f",secured[2][1]/1000))," [Gf/s] achieved by ",ApplyFunc(base_objid,Concat(secured[1],[tags])),"\n" );
           if (estimate[2][3]<=secured[2][1]) then
               over:=true;
           else
               Print("Best Hope is ",ApplyFunc(base_objid, Concat(estimate[1],[tags]))," with ",StringDouble("%.3f",estimate[2][1]/1000)," [Gf/s] (",StringDouble("%.3f",estimate[2][2]/1000)," < x < ",StringDouble("%.3f",estimate[2][3]/1000),")\n" );
               base_bench.exp.bases.benchTransforms:=[ApplyFunc(base_objid,Concat(estimate[1],[tags]))];
               base_bench.resumeAll();
               libgen_bench.runExhaust([ApplyFunc(libgen_objid,Concat(estimate[1]*2,[ApplyFunc(base_objid,Concat(estimate[1],[tags]))]))]);
               base_bench.exp.bases.benchTransforms:=[];
               base_bench.resumeAll(); #Fire Callbacks
           fi;
       od;
       return ApplyFunc(base_objid,Concat(secured[1],[tags]));
    end,
));


text:
#############################################################################
####
##
#W  anupqhead.gi               ANUPQ package                    Werner Nickel
#W                                                                Greg Gamble
##
##  `Head' file for the GAP interface to the ANU pq binary by Eamonn O'Brien.
##
#Y  Copyright (C) 2006  Lehrstuhl D fuer Mathematik,  RWTH Aachen,  Germany
##

#############################################################################
##
#V  ANUPQData . . record used by various functions of the ANUPQ package
##
##  The fields of ANUPQData are:
##
##    "binary"  . . the path of the pq binary
##    "tmpdir"  . . the path of the temporary directory for pq i/o files
##    "io"  . . . . list of data records for PqStart IO Streams
##    "outfile" . . the path of the pq output file
##    "SPimages"  . the path of the pq GAP_library file
##    "version" . . the version of the current pq binary
##
InstallValue( ANUPQData,
  rec( binary := Filename( DirectoriesPackagePrograms( "anupq" ), "pq" ),
       tmpdir := DirectoryTemporary(),
       ni := rec(), # record for non-interactive functions
       io := []     # list of records for PqStart IO Streams,
                    #  of which, there are initially none
       )
);
ANUPQData.outfile  := Filename( ANUPQData.tmpdir, "PQ_OUTPUT" );
ANUPQData.SPimages := Filename( ANUPQData.tmpdir, "GAP_library" );

# Fire up the pq binary to get its version
ANUPQData.version := "";
Process( DirectoryCurrent(), ANUPQData.binary, InputTextNone(),
         OutputTextString( ANUPQData.version, false ), [ "-v" ] );
ANUPQData.version :=
    ANUPQData.version{[PositionSublist( ANUPQData.version, "Version" ) + 8 ..
                       Length(ANUPQData.version) - 1] };

#############################################################################
##
#I  InfoClass
##
# Set the default level of InfoANUPQ
SetInfoLevel( InfoANUPQ, 1 );

#############################################################################
##
#V  ANUPQWarnOfOtherOptions . if true user is warned of non-ANUPQ-f'n options
##
ANUPQWarnOfOtherOptions := false;

#############################################################################
##
##  Ensure no zombie `pq' processes from interactive (`PqStart') sessions are
##  left lying around when user quits GAP.
##
InstallAtExit( PqQuitAll );

#E  anupqhead.gi . . . . . . . . . . . . . . . . . . . . . . . . .  ends here

max_stars_repo_name: gap-system/anupq

max_issues_repo_name: gap-system/anupq

max_forks_repo_name: gap-system/anupq


text:
--
-- An instance of this template must have a $Export section and the export_terminals option
--
-- Macros that may be redefined in an instance of this template
--
--     $eof_token
--     $additional_interfaces
--     $super_stream_class -- subclass com.ibm.lpg.LpgLexStream for getKind
--     $prs_stream_class -- use /.PrsStream./ if not subclassing
--
-- B E G I N N I N G   O F   T E M P L A T E   LexerTemplateD
--
%options programming_language=java,margin=4
%options table
%options action-block=("*.java", "/.", "./")
%options ParseTable=lpg.runtime.ParseTable
%options prefix=Char_

--
-- This template requires that the name of the EOF token be set
-- to EOF and that the prefix be "Char_" to be consistent with
-- KeywordTemplateD.
--
%Eof
    EOF
%End

--
-- This template also requires that the name of the parser EOF
-- Token to be exported be set to EOF_TOKEN
--
%Export
    EOF_TOKEN
%End

%Define
    --
    -- Macros that are be needed in an instance of this template
    --
    $eof_token /.$_EOF_TOKEN./

    $additional_interfaces /../
    $super_stream_class /.AbstractLexer./
    $prs_stream_class /.AbstractParser./
    $environment_class /.BasicEnvironment./
    $adapt_environment /.environment./

    $prs_stream /. // macro prs_stream is deprecated. Use function getPrsStream
                  getPrsStream()./
    $setSym1 /. // macro setSym1 is deprecated. Use function setResult
               lexParser.setSym1./
    $setResult /. // macro setResult is deprecated. Use function setResult
                 lexParser.setSym1./
    $getSym /. // macro getSym is deprecated. Use function getLastToken
              lexParser.getSym./
    $getToken /. // macro getToken is deprecated. Use function getToken
                lexParser.getToken./
    $getLeftSpan /. // macro getLeftSpan is deprecated. Use function getLeftSpan
                   lexParser.getFirstToken./
    $getRightSpan /. // macro getRightSpan is deprecated. Use function getRightSpan
                    lexParser.getLastToken./

    --
    -- Macros useful for specifying actions
    --
    $Header
    /.
                //
                // Rule $rule_number:  $rule_text
                //./

    $DefaultAction
    /. $Header
                case $rule_number: { ./

    $BeginAction /.$DefaultAction./

    $EndAction
    /.          break;
                }./

    $BeginJava
    /.$BeginAction
                $symbol_declarations./

    $EndJava /.$EndAction./

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

    $BeginActions
    /.
        public void ruleAction( int ruleNumber)
        {
            switch(ruleNumber)
            {./

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

	    public void ruleAction$rule_number(int ruleNumber)
	    {
	        switch (ruleNumber)
	        {./

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

%Globals
    /.import lpg.runtime.*;
    import org.eclipse.ocl.lpg.AbstractLexer;
    import org.eclipse.ocl.lpg.AbstractParser;
    ./
%End

%Headers
    /.
    @SuppressWarnings("nls")
    public class $action_type extends $super_stream_class implements $exp_type, $sym_type, RuleAction$additional_interfaces
    {
        private static ParseTable prs = new $prs_type();
        //
        // The Lexer contains an array of characters as the input stream to be parsed.
        // There are methods to retrieve and classify characters.
        // The lexparser "token" is implemented simply as the index of the next character in the array.
        // The Lexer extends the abstract class LpgLexStream with an implementation of the abstract
        // method getKind.  The template defines the Lexer class and the lexer() method.
        // A driver creates the action class, "Lexer", passing an Option object to the constructor.
        //
        protected $kw_lexer_class kwLexer;
        protected boolean printTokens;
        private $prs_stream_class parser;
        private LexParser lexParser = new LexParser(this, prs, this);

        private final $environment_class oclEnvironment;

        public $action_type($environment_class environment) {
            super($adapt_environment);
            oclEnvironment = environment;
        }

		public $action_class($environment_class environment, char[] chars) {
			this(environment, chars, "OCL", ECLIPSE_TAB_VALUE);
			kwLexer = new $kw_lexer_class(getInputChars(), $_IDENTIFIER);
		}

        public $action_type($environment_class environment, char[] input_chars, String filename, int tab)  {
            super($adapt_environment, input_chars, filename, tab);
            oclEnvironment = environment;
        }

		public $environment_class getOCLEnvironment() {
        	return oclEnvironment;
        }

        @Override
        public int [] getKeywordKinds() { return kwLexer.getKeywordKinds(); }

        public int getLeftSpan() { return lexParser.getFirstToken(); }
        public $prs_stream_class getParser() { return parser; }
        public int getRhsFirstTokenIndex(int i) { return lexParser.getFirstToken(i); }
        public int getRhsLastTokenIndex(int i) { return lexParser.getLastToken(i); }
        public int getRightSpan() { return lexParser.getLastToken(); }

        @Override
        public int getToken(int i) { return lexParser.getToken(i); }

        @Override
        public void initialize(char [] content, String filename)
        {
            super.initialize(content, filename);
            if (kwLexer == null)
                 kwLexer = new $kw_lexer_class(getInputChars(), $_IDENTIFIER);
            else
                 kwLexer.setInputChars(getInputChars());
        }

        @Override
        public String[] orderedExportedSymbols() { return $exp_type.orderedTerminalSymbols; }

	    @Override
	    public void setInputChars(char[] inputChars) {
			super.setInputChars(inputChars);
			kwLexer = new $kw_lexer_class(getInputChars(), $_IDENTIFIER);
		}

        @Override
        public void lexToTokens(Monitor monitor, $prs_stream_class parser)
        {
            if (getInputChars() == null)
                throw new NullPointerException("LexStream was not initialized");

            this.parser = parser;

            parser.makeToken(0, 0, 0); // Token list must start with a bad token

            lexParser.parseCharacters(monitor);  // Lex the input characters

            int i = getStreamIndex();
            parser.makeToken(i, i, $eof_token); // and end with the end of file token
            parser.setStreamLength(parser.getSize());

            return;
        }
    ./
%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

_AVXINTRIN := () ->  ["<immintrin.h>"];
_MM_MALLOC := () -> When(not LocalConfig.osinfo.isDarwin(), ["<include/mm_malloc.h>"],[]);

Class(AVX_Intel, SIMD_ISA, rec(
    unparser := AVXUnparser,
    arch := "Intel_AVX",
    info := "Intel AVX architecture definition base",
    declareConstants := false,

    file := "avx",
    #arrayDataModifier := "static __declspec(align(32))",
    #arrayBufModifier  := "static __declspec(align(32))",
	alignmentBytes := 32,

    autolib := rec(
        includes := () -> _NMMINTRIN() :: _SMMINTRIN() :: _TMMINTRIN() :: _PMMINTRIN()
                       :: _EMMINTRIN() :: _XMMINTRIN() :: _AVXINTRIN(),
        timerIncludes := () -> ["<include/sp_rdtsc.h>"],
        cFlags := "-mavx",
    ),

    simpIndicesInside := AVX_LDST.getList(),
));

#==============================================================================================
# 4-way double precision real

Class(AVX_4x64f, AVX_Intel, rec(
    info     := "AVX 4 x 64-bit double",
    countrec := rec(
        ops := [
            [ add, sub, addsub_4x64f, chslo_2x64f, chshi_2x64f, addsub_2x64f, hadd_2x64f ],
	    [ mul ],
	    [ fmaddsub_4x64f],
            [ vunpacklo_4x64f,  vunpackhi_4x64f,  vshuffle_4x64f,  vperm2_4x64f,    vpermf128_4x64f,
              vuunpacklo_4x64f, vuunpackhi_4x64f, vushuffle_4x64f, vuperm2_4x64f,   vupermf128_4x64f,
              vunpacklo_2x64f,  vunpackhi_2x64f,  vshuffle_2x64f,  vushuffle_2x64f, vperm_4x64f,
              vinsert_2l_4x64f, vextract_2l_4x64f, vblend_4x64f],
            [ vloadmask_4x64f,  vstoremask_4x64f, vload1sd_2x64f,  vload_1l_2x64f,  vload_1h_2x64f,
	      vloadu_2x64f,     vstore_1l_2x64f,  vstore_1h_2x64f, vstoreu_2x64f ],
            [ deref, vloadu_4x64f, vstoreu_4x64f ],
            Value      # Value without [] is a keyword in countOps !!
        ],
        printstrings := ["[adds]", "[mults]", "[fmas]", "[vperms]", "[svldst]", "[vldst]", "[vval]"],
        arithcost := (self, opcount) >> opcount[1]+opcount[2]+2*opcount[3]
    ),

    includes     := () -> ["<include/omega64.h>"] :: _AVXINTRIN(),
    active       := true,
    isFixedPoint := false,
    isFloat      := true,

    v       := 4,
    t       := TVect(T_Real(64), 4),
    ctype   := "double",
    instr   := [ vunpacklo_4x64f, vunpackhi_4x64f, vshuffle_4x64f, vperm2_4x64f, vpermf128_4x64f, vperm_4x64f, vblend_4x64f ],
               #vuunpacklo_4x64f, vuunpackhi_4x64f, vushuffle_4x64f, vuperm2_4x64f, vupermf128_4x64f
    bits    := 64,
    splopts := rec(precision := "double"),

    freshU  := self >> var.fresh_t("U", self.t),

    dupload := (y, x) -> assign(y, vdup(x, 4)),
    mul_cx := (self, opts) >>
        ((y, x, c) -> let( u1 := self.freshU(), u2 := self.freshU(), u3 := self.freshU(),
            decl([u1, u2, u3], chain(
                assign(u1, mul(x, vunpacklo_4x64f(c, c))),
                assign(u2, vshuffle_4x64f(x, x, [2,1,2,1])),
                assign(u3, mul(u2, vunpackhi_4x64f(c, c))),
                assign(y, addsub_4x64f(u1, u3))))
        )),

    mul_cx_conj := (self, opts) >>
        ((y,x,c) -> let(
		u1 := self.freshU(), u2 := self.freshU(), v := self.freshU(),
                decl([u1, u2, v], chain(
                    assign(u1,             mul(x, vunpackhi_4x64f(c, c))),
                    assign(u2, mul(vshuffle_4x64f(x, x, [2,1,2,1]), vunpacklo_4x64f(c, c))),
                    assign(v, addsub_4x64f(u2, u1)),
                    assign(y, vshuffle_4x64f(v, v, [2,1,2,1]))
                )))),

    svload_init := (vt) -> [
        # load using subvectors of length 1
	#
	[
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 2)),
            decl([u1], chain(
                assign(u1, vload1sd_2x64f(x[1].toPtr(vt.t))),
                assign(y, vinsert_2l_4x64f(vt.zero(), u1, [0]))
            ))),
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 2)), u2 := var.fresh_t("U", TVect(vt.t, 2)),
            decl([u1, u2], chain(
                assign(u1, vload1sd_2x64f(x[1].toPtr(vt.t))),
                assign(u2, vload_1h_2x64f(u1, x[2].toPtr(vt.t))),
                assign(y, vinsert_2l_4x64f(vt.zero(), u2, [0]))
            ))),
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 2)), u2 := var.fresh_t("U", TVect(vt.t, 2)),
                        u3 := var.fresh_t("U", TVect(vt.t, 2)), v1 := var.fresh_t("U", vt),
            decl([u1, u2, u3, v1], chain(
                assign(u1, vload1sd_2x64f(x[1].toPtr(vt.t))),
                assign(u2, vload_1h_2x64f(u1, x[2].toPtr(vt.t))),
                assign(u3, vload1sd_2x64f(x[3].toPtr(vt.t))),
                assign(v1, vinsert_2l_4x64f(vt.zero(), u2, [0])),
                assign(y, vinsert_2l_4x64f(v1, u3, [1]))
            ))),
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 2)), u2 := var.fresh_t("U", TVect(vt.t, 2)),
                        u3 := var.fresh_t("U", TVect(vt.t, 2)), u4 := var.fresh_t("U", TVect(vt.t, 2)),
                        v1 := var.fresh_t("U", vt),
            decl([u1, u2, u3, v1], chain(
                assign(u1, vload1sd_2x64f(x[1].toPtr(vt.t))),
                assign(u2, vload_1h_2x64f(u1, x[2].toPtr(vt.t))),
                assign(u3, vload1sd_2x64f(x[3].toPtr(vt.t))),
                assign(u4, vload_1h_2x64f(u3, x[4].toPtr(vt.t))),
                assign(v1, vinsert_2l_4x64f(vt.zero(), u2, [0])),
                assign(y, vinsert_2l_4x64f(v1, u4, [1]))
            )))
        ],
        # load using subvectors of length 2
	#
        [
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 2)),
            decl([u1], chain(
                assign(u1, vloadu_2x64f(x[1].toPtr(vt.t))),
                assign(y, vinsert_2l_4x64f(vt.zero(), u1, [0]))
            ))),
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 2)), u2 := var.fresh_t("U", TVect(vt.t, 2)),
                        v1 := var.fresh_t("U", vt),
            decl([u1, u2, v1], chain(
                assign(u1, vloadu_2x64f(x[1].toPtr(vt.t))),
                assign(u2, vloadu_2x64f(x[2].toPtr(vt.t))),
                assign(v1, vinsert_2l_4x64f(vt.zero(), u1, [0])),
                assign(y, vinsert_2l_4x64f(v1, u2, [1]))
            )))
        ]
    ],
    svstore_init := (vt) -> [
        # store using subvectors of length 1
        [
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 2)),
            decl([u1], chain(
                assign(u1, vextract_2l_4x64f(x, [0])),
                vstore_1l_2x64f(y[1].toPtr(vt.t), u1)
            ))),
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 2)),
            decl([u1], chain(
                assign(u1, vextract_2l_4x64f(x, [0])),
                vstore_1l_2x64f(y[1].toPtr(vt.t), u1),
                vstore_1h_2x64f(y[2].toPtr(vt.t), u1)
            ))),
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 2)), u2 := var.fresh_t("U", TVect(vt.t, 2)),
            decl([u1, u2], chain(
                assign(u1, vextract_2l_4x64f(x, [0])),
                vstore_1l_2x64f(y[1].toPtr(vt.t), u1),
                vstore_1h_2x64f(y[2].toPtr(vt.t), u1),
                assign(u2, vextract_2l_4x64f(x, [1])),
                vstore_1l_2x64f(y[3].toPtr(vt.t), u2)
            ))),
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 2)), u2 := var.fresh_t("U", TVect(vt.t, 2)),
            decl([u1, u2], chain(
                assign(u1, vextract_2l_4x64f(x, [0])),
                vstore_1l_2x64f(y[1].toPtr(vt.t), u1),
                vstore_1h_2x64f(y[2].toPtr(vt.t), u1),
                assign(u2, vextract_2l_4x64f(x, [1])),
                vstore_1l_2x64f(y[3].toPtr(vt.t), u2),
                vstore_1h_2x64f(y[4].toPtr(vt.t), u2)
            )))
        ],
        # store using subvectors of length 2
        [
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 2)),
            decl([u1], chain(
                assign(u1, vextract_2l_4x64f(x, [0])),
                vstoreu_2x64f(y[1].toPtr(vt.t), u1)
            ))),
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 2)), u2 := var.fresh_t("U", TVect(vt.t, 2)),
            decl([u1, u2], chain(
                assign(u1, vextract_2l_4x64f(x, [0])),
                vstoreu_2x64f(y[1].toPtr(vt.t), u1),
                assign(u2, vextract_2l_4x64f(x, [1])),
                vstoreu_2x64f(y[2].toPtr(vt.t), u2)
            )))
        ]],

    svload  := ~.svload_init(~.t),
    svstore := ~.svstore_init(~.t),

    loadc := (self, sv, opts) >> Cond(
        sv = 4,
            ((y,x) -> assign(y, vloadu_4x64f(x.toPtr(self.t.t)))),
        sv = 2,
            ((y,x) -> assign(y, vinsert_2l_4x64f(self.t.zero(), vloadu_2x64f(x.toPtr(self.t.t)), [0]))),
        # else
	Checked( sv=1,
	    ((y,x) -> assign(y, vinsert_2l_4x64f(self.t.zero(), vload1sd_2x64f(x.toPtr(self.t.t)), [0]))))),

    loadc_align := (self, sv, align, opts) >> Cond(
        align = 0 and sv = 8,
            ((y,x,addr) -> assign(y, deref(nth(x,addr).toPtr(self.t)))),
        align mod 2 = 0 and sv = 2,
            ((y,x,addr) -> assign(y, vinsert_2l_4x64f(self.t.zero(), deref(nth(x,addr).toPtr(TVect(self.t.t, 2))), [0]))),
        sv = 4,
            ((y,x,addr) -> assign(y, vloadu_4x64f(nth(x, addr).toPtr(self.t.t)))),
        sv = 2,
            ((y,x,addr) -> assign(y, vinsert_2l_4x64f(self.t.zero(), vloadu_2x64f(nth(x, addr).toPtr(self.t.t)), [0]))),
        # else
        Checked( sv=1,
	    ((y,x,addr) -> assign(y, vinsert_2l_4x64f(self.t.zero(), vload1sd_2x64f(nth(x, addr).toPtr(self.t.t)), [0]))))),

    storec_init := (vt) -> [
        ((y,x) -> vstore_1l_2x64f(y.toPtr(vt.t), vextract_2l_4x64f(x, [0]))),
        ((y,x) -> vstoreu_2x64f(y.toPtr(vt.t), vextract_2l_4x64f(x, [0]))),
        ((y,x) -> vstoremask_4x64f(y.toPtr(vt.t), x, vstoremask_4x64f.mask(3))),
        ((y,x) -> vstoreu_4x64f(y.toPtr(vt.t), x))
    ],

    storec := ~.storec_init(~.t),

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

    bin_shrev := (self, y, x, opts) >> let(
        u1 := self.freshU(), u2 := self.freshU(),
        decl( [u1, u2], chain(
            assign(u1, vupermf128_4x64f(x[2], [2,1])),
            assign(u2, vblend_4x64f(x[1], x[2], [1,1,2,2])),
            assign(y,  vblend_4x64f(u1, u2, [2,1,2,1]))))),
));

#==============================================================================================
Class(AVX_8x32f, AVX_Intel, rec(
    info     := "AVX 8 x 32-bit float",
    countrec := rec(
        ops := [
            [ add, sub, addsub_8x32f, addsub_4x32f, hadd_4x32f ],
	    [ mul ],
	    [ fmaddsub_8x32f ],
            [ vunpacklo_8x32f,  vunpackhi_8x32f,   vpermf128_8x32f,  vperm2_8x32f,  vshuffle_8x32f,
              vinsert_4l_8x32f, vextract_4l_8x32f, vperm_8x32f,      vpermv_8x32f,
              vunpacklo_4x32f,  vunpackhi_4x32f,   vshuffle_4x32f,   vushuffle_4x32f,
              vuunpacklo_8x32f, vuunpackhi_8x32f,  vupermf128_8x32f, vushuffle_8x32f,
              vldup_8x32f, vhdup_8x32f, vblend_8x32f ],
            [ vloadmask_8x32f, vstoremask_8x32f,
              vload1_4x32f,    vload_2l_4x32f,  vload_2h_4x32f,  vloadu_4x32f,  vloadu2_4x32f,
              vstore1_4x32f,   vstore_2l_4x32f, vstore_2h_4x32f, vstoreu_4x32f, vstoreu2_4x32f,
              vstoremsk_4x32f ],
            [ deref, vloadu_8x32f, vstoreu_8x32f ],
            Value   # Value without [] is a keyword in countOps !!
        ],
        printstrings := ["[adds]", "[mults]", "[fmas]", "[vperms]", "[svldst]", "[vldst]", "[vval]"],
        arithcost := (self, opcount) >> opcount[1]+opcount[2]+2*opcount[3]
    ),

    includes     := () -> ["<include/omega32.h>"] :: _AVXINTRIN(),
    active       := true,
    isFixedPoint := false,
    isFloat      := true,

    v       := 8,
    t       := TVect(T_Real(32), 8),
    ctype   := "float",
    instr   := [ vunpacklo_8x32f,  vunpackhi_8x32f,  vpermf128_8x32f,  vshuffle_8x32f,  vperm_8x32f,
                 vuunpacklo_8x32f, vuunpackhi_8x32f, vupermf128_8x32f, vushuffle_8x32f, vblend_8x32f,
                 vcxtr_8x32f],

    bits    := 32,
    splopts := rec(precision := "single"),

    freshU  := self >> var.fresh_t("U", self.t),

    dupload := (y, x) -> assign(y, vdup(x, 8)),

    mul_cx := (self, opts) >>
        ((y, x, c) -> let( u1 := self.freshU(), u2 := self.freshU(), u3 := self.freshU(),
            decl([u1, u2, u3], chain(
                assign(u1, mul(x, vldup_8x32f(c))),
                assign(u2, vshuffle_8x32f(x, x, [2,1,4,3])),
                assign(u3, mul(u2, vhdup_8x32f(c))),
                assign(y,  addsub_8x32f(u1, u3))))
        )),

    mul_cx_conj := (self, opts) >>
        ((y,x,c) -> let( u1 := self.freshU(), u2 := self.freshU(), v := self.freshU(),
            decl([u1, u2, v], chain(
                assign(u1, mul(x, vhdup_8x32f(c))),
                assign(u2, mul(vshuffle_8x32f(x, x, [2,1,4,3]), vldup_8x32f(c))),
                assign(v, addsub_8x32f(u2, u1)),
                assign(y, vshuffle_8x32f(v, v, [2,1,4,3]))))
        )),

    svload_init := (vt) -> [
        # Load using subvectors of length 1.
	#
	[
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 4)),
            decl([u1], chain(
                assign(u1, vload1_4x32f(x[1].toPtr(vt.t))),
                assign(y,  vinsert_4l_8x32f(vt.zero(), u1, [0]))
            ))),
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 4)), u2 := var.fresh_t("U", TVect(vt.t, 4)),
            decl([u1, u2], chain(
                assign(u1, vload1_4x32f(x[1].toPtr(vt.t))),
                assign(u2, vinsert_4x32f(u1, deref(x[2].toPtr(T_Int(32))), 2)),
                assign(y,  vinsert_4l_8x32f(vt.zero(), u2, [0]))
            ))),
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 4)), u2 := var.fresh_t("U", TVect(vt.t, 4)),
                     u3 := var.fresh_t("U", TVect(vt.t, 4)),
            decl([u1, u2, u3], chain(
                assign(u1, vload1_4x32f(x[1].toPtr(vt.t))),
                assign(u2, vinsert_4x32f(u1, deref(x[2].toPtr(T_Int(32))), 2)),
                assign(u3, vinsert_4x32f(u2, deref(x[3].toPtr(T_Int(32))), 3)),
                assign(y,  vinsert_4l_8x32f(vt.zero(), u3, [0]))
            ))),
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 4)), u2 := var.fresh_t("U", TVect(vt.t, 4)),
                     u3 := var.fresh_t("U", TVect(vt.t, 4)), u4 := var.fresh_t("U", TVect(vt.t, 4)),
            decl([u1, u2, u3, u4], chain(
                assign(u1, vload1_4x32f(x[1].toPtr(vt.t))),
                assign(u2, vinsert_4x32f(u1, deref(x[2].toPtr(T_Int(32))), 2)),
                assign(u3, vinsert_4x32f(u2, deref(x[3].toPtr(T_Int(32))), 3)),
                assign(u4, vinsert_4x32f(u3, deref(x[4].toPtr(T_Int(32))), 4)),
                assign(y,  vinsert_4l_8x32f(vt.zero(), u4, [0]))
            ))),
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 4)), u2 := var.fresh_t("U", TVect(vt.t, 4)),
                     u3 := var.fresh_t("U", TVect(vt.t, 4)), u4 := var.fresh_t("U", TVect(vt.t, 4)),
                     u5 := var.fresh_t("U", TVect(vt.t, 4)),
                     v1 := var.fresh_t("U", vt),
            decl([u1, u2, u3, u4, u5, v1], chain(
                assign(u1, vload1_4x32f(x[1].toPtr(vt.t))),
                assign(u2, vinsert_4x32f(u1, deref(x[2].toPtr(T_Int(32))), 2)),
                assign(u3, vinsert_4x32f(u2, deref(x[3].toPtr(T_Int(32))), 3)),
                assign(u4, vinsert_4x32f(u3, deref(x[4].toPtr(T_Int(32))), 4)),
                assign(v1, vinsert_4l_8x32f(vt.zero(), u4, [0])),
                assign(u5, vload1_4x32f(x[5].toPtr(vt.t))),
                assign(y,  vinsert_4l_8x32f(v1, u5, [1]))
            ))),

        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 4)), u2 := var.fresh_t("U", TVect(vt.t, 4)),
                     u3 := var.fresh_t("U", TVect(vt.t, 4)), u4 := var.fresh_t("U", TVect(vt.t, 4)),
                     u5 := var.fresh_t("U", TVect(vt.t, 4)), u6 := var.fresh_t("U", TVect(vt.t, 4)),
                     v1 := var.fresh_t("U", vt),
            decl([u1, u2, u3, u4, u5, u6, v1], chain(
                assign(u1, vload1_4x32f(x[1].toPtr(vt.t))),
                assign(u2, vinsert_4x32f(u1, deref(x[2].toPtr(T_Int(32))), 2)),
                assign(u3, vinsert_4x32f(u2, deref(x[3].toPtr(T_Int(32))), 3)),
                assign(u4, vinsert_4x32f(u3, deref(x[4].toPtr(T_Int(32))), 4)),
                assign(v1, vinsert_4l_8x32f(vt.zero(), u4, [0])),
                assign(u5, vload1_4x32f(x[5].toPtr(vt.t))),
                assign(u6, vinsert_4x32f(u5, deref(x[6].toPtr(T_Int(32))), 2)),
                assign(y,  vinsert_4l_8x32f(v1, u6, [1]))
            ))),
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 4)), u2 := var.fresh_t("U", TVect(vt.t, 4)),
                     u3 := var.fresh_t("U", TVect(vt.t, 4)), u4 := var.fresh_t("U", TVect(vt.t, 4)),
                     u5 := var.fresh_t("U", TVect(vt.t, 4)), u6 := var.fresh_t("U", TVect(vt.t, 4)),
                     u7 := var.fresh_t("U", TVect(vt.t, 4)),
                     v1 := var.fresh_t("U", vt),
            decl([u1, u2, u3, u4, u5, u6, u7, v1], chain(
                assign(u1, vload1_4x32f(x[1].toPtr(vt.t))),
                assign(u2, vinsert_4x32f(u1, deref(x[2].toPtr(T_Int(32))), 2)),
                assign(u3, vinsert_4x32f(u2, deref(x[3].toPtr(T_Int(32))), 3)),
                assign(u4, vinsert_4x32f(u3, deref(x[4].toPtr(T_Int(32))), 4)),
                assign(v1, vinsert_4l_8x32f(vt.zero(), u4, [0])),
                assign(u5, vload1_4x32f(x[5].toPtr(vt.t))),
                assign(u6, vinsert_4x32f(u5, deref(x[6].toPtr(T_Int(32))), 2)),
                assign(u7, vinsert_4x32f(u6, deref(x[7].toPtr(T_Int(32))), 3)),
                assign(y,  vinsert_4l_8x32f(v1, u7, [1]))
            ))),
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 4)), u2 := var.fresh_t("U", TVect(vt.t, 4)),
                     u3 := var.fresh_t("U", TVect(vt.t, 4)), u4 := var.fresh_t("U", TVect(vt.t, 4)),
                     u5 := var.fresh_t("U", TVect(vt.t, 4)), u6 := var.fresh_t("U", TVect(vt.t, 4)),
                     u7 := var.fresh_t("U", TVect(vt.t, 4)), u8 := var.fresh_t("U", TVect(vt.t, 4)),
                     v1 := var.fresh_t("U", vt),
            decl([u1, u2, u3, u4, u5, u6, u7, u8, v1], chain(
                assign(u1, vload1_4x32f(x[1].toPtr(vt.t))),
                assign(u2, vinsert_4x32f(u1, deref(x[2].toPtr(T_Int(32))), 2)),
                assign(u3, vinsert_4x32f(u2, deref(x[3].toPtr(T_Int(32))), 3)),
                assign(u4, vinsert_4x32f(u3, deref(x[4].toPtr(T_Int(32))), 4)),
                assign(v1, vinsert_4l_8x32f(vt.zero(), u4, [0])),
                assign(u5, vload1_4x32f(x[5].toPtr(vt.t))),
                assign(u6, vinsert_4x32f(u5, deref(x[6].toPtr(T_Int(32))), 2)),
                assign(u7, vinsert_4x32f(u6, deref(x[7].toPtr(T_Int(32))), 3)),
                assign(u8, vinsert_4x32f(u7, deref(x[8].toPtr(T_Int(32))), 4)),
                assign(y,  vinsert_4l_8x32f(v1, u8, [1]))
            )))
        ],
        # Load using subvectors of length 2
	#
        [
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 4)),
            decl([u1], chain(
                assign(u1, vload_2l_4x32f(vzero_4x32f(), x[1].toPtr(TVect(vt.t, 2)))),
                assign(y,  vinsert_4l_8x32f(vt.zero(), u1, [0]))
            ))),
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 4)), u2 := var.fresh_t("U", TVect(vt.t, 4)),
            decl([u1, u2], chain(
                assign(u1, vload_2l_4x32f(vzero_4x32f(), x[1].toPtr(TVect(vt.t, 2)))),
                assign(u2, vload_2h_4x32f(u1, x[2].toPtr(TVect(vt.t, 2)))),
                assign(y,  vinsert_4l_8x32f(vt.zero(), u2, [0]))
            ))),
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 4)), u2 := var.fresh_t("U", TVect(vt.t, 4)),
                     u3 := var.fresh_t("U", TVect(vt.t, 4)),
                     z1 := var.fresh_t("Z", TVect(vt.t, 4)),
                     v1 := var.fresh_t("U", vt),
            decl([u1, u2, u3, v1], chain(
                assign(z1, vzero_4x32f()),
                assign(u1, vload_2l_4x32f(z1, x[1].toPtr(TVect(vt.t, 2)))),
                assign(u2, vload_2h_4x32f(u1, x[2].toPtr(TVect(vt.t, 2)))),
                assign(v1, vinsert_4l_8x32f(vt.zero(), u2, [0])),
                assign(u3, vload_2l_4x32f(z1, x[3].toPtr(TVect(vt.t, 2)))),
                assign(y,  vinsert_4l_8x32f(v1, u3, [1]))
            ))),
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 4)), u2 := var.fresh_t("U", TVect(vt.t, 4)),
                     u3 := var.fresh_t("U", TVect(vt.t, 4)), u4 := var.fresh_t("U", TVect(vt.t, 4)),
                     z1 := var.fresh_t("Z", TVect(vt.t, 4)),
                     v1 := var.fresh_t("U", vt),
            decl([u1, u2, u3, u4, v1], chain(
                assign(z1, vzero_4x32f()),
                assign(u1, vload_2l_4x32f(z1, x[1].toPtr(TVect(vt.t, 2)))),
                assign(u2, vload_2h_4x32f(u1, x[2].toPtr(TVect(vt.t, 2)))),
                assign(v1, vinsert_4l_8x32f(vt.zero(), u2, [0])),
                assign(u3, vload_2l_4x32f(z1, x[3].toPtr(TVect(vt.t, 2)))),
                assign(u4, vload_2h_4x32f(u3, x[4].toPtr(TVect(vt.t, 2)))),
                assign(y,  vinsert_4l_8x32f(v1, u4, [1]))
            ))),
        ],
        # Load using subvectors of length 4.
	#
        [
        (y,x,opts) ->
            assign(y, vinsert_4l_8x32f(vt.zero(), vloadu_4x32f(x[1].toPtr(vt.t)), [0])),

        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 4)),
            decl([u1], chain(
                assign(u1, vinsert_4l_8x32f(vt.zero(), vloadu_4x32f(x[1].toPtr(vt.t)), [0])),
                assign(y,  vinsert_4l_8x32f(u1,        vloadu_4x32f(x[2].toPtr(vt.t)), [1]))
            )))
        ]
    ],

    svstore_init := (vt) -> [
        # Store using subvectors of length 1.
	#
        [
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 4)),
            decl([u1], chain(
                assign(u1, vextract_4l_8x32f(x, [0])),
                vstore1_4x32f(y[1].toPtr(vt.t), u1)
            ))),
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 4)),
            decl([u1], chain(
                assign(u1, vextract_4l_8x32f(x, [0])),
                vstore1_4x32f(y[1].toPtr(vt.t), u1),
                vextract_4x32f(y[2].toPtr(T_Int(32)), u1, 2)
            ))),
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 4)),
            decl([u1], chain(
                assign(u1, vextract_4l_8x32f(x, [0])),
                vstore1_4x32f(y[1].toPtr(vt.t), u1),
                vextract_4x32f(y[2].toPtr(T_Int(32)), u1, 2),
                vextract_4x32f(y[3].toPtr(T_Int(32)), u1, 3)
            ))),
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 4)),
            decl([u1], chain(
                assign(u1, vextract_4l_8x32f(x, [0])),
                vstore1_4x32f(y[1].toPtr(vt.t), u1),
                vextract_4x32f(y[2].toPtr(T_Int(32)), u1, 2),
                vextract_4x32f(y[3].toPtr(T_Int(32)), u1, 3),
                vextract_4x32f(y[4].toPtr(T_Int(32)), u1, 4)
            ))),

        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 4)), u2 := var.fresh_t("U", TVect(vt.t, 4)),
            decl([u1, u2], chain(
                assign(u1, vextract_4l_8x32f(x, [0])),
                vstore1_4x32f(y[1].toPtr(vt.t), u1),
                vextract_4x32f(y[2].toPtr(T_Int(32)), u1, 2),
                vextract_4x32f(y[3].toPtr(T_Int(32)), u1, 3),
                vextract_4x32f(y[4].toPtr(T_Int(32)), u1, 4),
                assign(u2, vextract_4l_8x32f(x, [1])),
                vstore1_4x32f(y[5].toPtr(vt.t), u2)
            ))),
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 4)), u2 := var.fresh_t("U", TVect(vt.t, 4)),
            decl([u1, u2], chain(
                assign(u1, vextract_4l_8x32f(x, [0])),
                vstore1_4x32f(y[1].toPtr(vt.t), u1),
                vextract_4x32f(y[2].toPtr(T_Int(32)), u1, 2),
                vextract_4x32f(y[3].toPtr(T_Int(32)), u1, 3),
                vextract_4x32f(y[4].toPtr(T_Int(32)), u1, 4),
                assign(u2, vextract_4l_8x32f(x, [1])),
                vstore1_4x32f(y[5].toPtr(vt.t), u2),
                vextract_4x32f(y[6].toPtr(T_Int(32)), u2, 2)
            ))),
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 4)), u2 := var.fresh_t("U", TVect(vt.t, 4)),
            decl([u1, u2], chain(
                assign(u1, vextract_4l_8x32f(x, [0])),
                vstore1_4x32f(y[1].toPtr(vt.t), u1),
                vextract_4x32f(y[2].toPtr(T_Int(32)), u1, 2),
                vextract_4x32f(y[3].toPtr(T_Int(32)), u1, 3),
                vextract_4x32f(y[4].toPtr(T_Int(32)), u1, 4),
                assign(u2, vextract_4l_8x32f(x, [1])),
                vstore1_4x32f(y[5].toPtr(vt.t), u2),
                vextract_4x32f(y[6].toPtr(T_Int(32)), u2, 2),
                vextract_4x32f(y[7].toPtr(T_Int(32)), u2, 3)
            ))),
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 4)), u2 := var.fresh_t("U", TVect(vt.t, 4)),
            decl([u1, u2], chain(
                assign(u1, vextract_4l_8x32f(x, [0])),
                vstore1_4x32f(y[1].toPtr(vt.t), u1),
                vextract_4x32f(y[2].toPtr(T_Int(32)), u1, 2),
                vextract_4x32f(y[3].toPtr(T_Int(32)), u1, 3),
                vextract_4x32f(y[4].toPtr(T_Int(32)), u1, 4),
                assign(u2, vextract_4l_8x32f(x, [1])),
                vstore1_4x32f(y[5].toPtr(vt.t), u2),
                vextract_4x32f(y[6].toPtr(T_Int(32)), u2, 2),
                vextract_4x32f(y[7].toPtr(T_Int(32)), u2, 3),
                vextract_4x32f(y[8].toPtr(T_Int(32)), u2, 4)
            )))
        ],
        # Store using subvectors of length 2.
	#
        [
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 4)),
            decl([u1], chain(
                assign(u1, vextract_4l_8x32f(x, [0])),
                vstore_2l_4x32f(y[1].toPtr(TVect(vt.t, 2)), u1)
            ))),
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 4)),
            decl([u1], chain(
                assign(u1, vextract_4l_8x32f(x, [0])),
                vstore_2l_4x32f(y[1].toPtr(TVect(vt.t, 2)), u1),
                vstore_2h_4x32f(y[2].toPtr(TVect(vt.t, 2)), u1)
            ))),
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 4)), u2 := var.fresh_t("U", TVect(vt.t, 4)),
            decl([u1, u2], chain(
                assign(u1, vextract_4l_8x32f(x, [0])),
                vstore_2l_4x32f(y[1].toPtr(TVect(vt.t, 2)), u1),
                vstore_2h_4x32f(y[2].toPtr(TVect(vt.t, 2)), u1),
                assign(u2, vextract_4l_8x32f(x, [1])),
                vstore_2l_4x32f(y[3].toPtr(TVect(vt.t, 2)), u2)
            ))),
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 4)), u2 := var.fresh_t("U", TVect(vt.t, 4)),
            decl([u1, u2], chain(
                assign(u1, vextract_4l_8x32f(x, [0])),
                vstore_2l_4x32f(y[1].toPtr(TVect(vt.t, 2)), u1),
                vstore_2h_4x32f(y[2].toPtr(TVect(vt.t, 2)), u1),
                assign(u2, vextract_4l_8x32f(x, [1])),
                vstore_2l_4x32f(y[3].toPtr(TVect(vt.t, 2)), u2),
                vstore_2h_4x32f(y[4].toPtr(TVect(vt.t, 2)), u2)
            )))
        ],
        # Store using subvectors of length 4.
	#
        [
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 4)),
            decl([u1], chain(
                assign(u1, vextract_4l_8x32f(x, [0])),
                vstoreu_4x32f(y[1].toPtr(vt.t), u1)
            ))),
        (y,x,opts) -> let(u1 := var.fresh_t("U", TVect(vt.t, 4)), u2 := var.fresh_t("U", TVect(vt.t, 4)),
            decl([u1, u2], chain(
                assign(u1, vextract_4l_8x32f(x, [0])),
                vstoreu_4x32f(y[1].toPtr(vt.t), u1),
                assign(u2, vextract_4l_8x32f(x, [1])),
                vstoreu_4x32f(y[2].toPtr(vt.t), u2)
            )))
        ]
    ],

    svload  := ~.svload_init(~.t),
    svstore := ~.svstore_init(~.t),

    loadc := (self, sv, opts) >> Cond(
        sv = 8,
            ((y,x) -> assign(y, vloadu_8x32f(x.toPtr(self.t.t)))),
        sv = 4,
            ((y,x) -> assign(y, vinsert_4l_8x32f(self.t.zero(), vloadu_4x32f(x.toPtr(self.t.t)), [0]))),
        sv = 2,
            ((y,x) -> assign(y, vinsert_4l_8x32f(self.t.zero(), vload_2l_4x32f(vzero_4x32f(), x.toPtr(TVect(self.t.t, 2))), [0]))),
        # else
            ((y,x) -> assign(y, vinsert_4l_8x32f(self.t.zero(), vload1_4x32f(x.toPtr(self.t.t)), [0])))),

    loadc_align := (self, sv, align, opts) >> Cond(
        align = 0 and sv = 8,
            ((y,x,addr) -> assign(y, deref(nth(x,addr).toPtr(self.t)))),
        align mod 4 = 0 and sv = 4,
            ((y,x,addr) -> assign(y, vinsert_4l_8x32f(self.t.zero(), deref(nth(x,addr).toPtr(TVect(self.t.t, 4))), [0]))),
        sv = 8,
            ((y,x,addr) -> assign(y, vloadu_8x32f(nth(x,addr).toPtr(self.t.t)))),
        sv = 4,
            ((y,x,addr) -> assign(y, vinsert_4l_8x32f(self.t.zero(), vloadu_4x32f(nth(x,addr).toPtr(self.t.t)), [0]))),
        sv = 2,
            ((y,x,addr) -> assign(y, vinsert_4l_8x32f(self.t.zero(), vload_2l_4x32f(vzero_4x32f(), nth(x,addr).toPtr(TVect(self.t.t, 2))), [0]))),
        # else
            ((y,x,addr) -> assign(y, vinsert_4l_8x32f(self.t.zero(), vload1_4x32f(nth(x,addr).toPtr(self.t.t)), [0])))),

    storec_init := (vt) -> [
        ((y,x) -> vstore1_4x32f(y.toPtr(vt.t), vextract_4l_8x32f(x, [0]))),
        ((y,x) -> vstore_2l_4x32f(y.toPtr(TVect(vt.t, 2)), vextract_4l_8x32f(x, [0]))),
        ((y,x) -> vstoremask_8x32f(y.toPtr(vt.t), x, vstoremask_8x32f.mask(3))),
        ((y,x) -> vstoreu_4x32f(y.toPtr(vt.t), vextract_4l_8x32f(x, [0]))),
        ((y,x) -> vstoremask_8x32f(y.toPtr(vt.t), x, vstoremask_8x32f.mask(5))),
        ((y,x) -> vstoremask_8x32f(y.toPtr(vt.t), x, vstoremask_8x32f.mask(6))),
        ((y,x) -> vstoremask_8x32f(y.toPtr(vt.t), x, vstoremask_8x32f.mask(7))),
        ((y,x) -> vstoreu_8x32f(y.toPtr(vt.t), x))
    ],

    storec  := ~.storec_init(~.t),

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

    # VJxI implements VJxI(m,v) which is J(m) /tensor I(v/m)
    VJxI := (m, y, x, opts) -> Cond( m=1, assign(y, x),
                                     m=2, assign(y, vupermf128_8x32f(x, [2,1])),
                                     m=4, assign(y, vperm_8x32f(vupermf128_8x32f(x, [2,1]), [3,4,1,2])),
                                     m=8, assign(y, vperm_8x32f(vupermf128_8x32f(x, [2,1]), [4,3,2,1])),
                                     Error("unexpected m")),

    bin_shrev := (self, y, x, opts) >> let(
        u1 := self.freshU(), u2 := self.freshU(), u3 := self.freshU(),
        decl( [u1, u2, u3], chain(
            assign(u1, vupermf128_8x32f(x[2], [2,1])),                # [a,b,c,d,e,f,g,h] -> [e,f,g,h,a,b,c,d]
            assign(u2, vblend_8x32f(x[1], x[2], [1,1,1,1,2,2,2,2])),  # -> [x11,x12,x13,x14,e,f,g,h]
            assign(u3, vushuffle_8x32f(u1, [1,4,3,2])),               # [e,f,g,h,a,b,c,d] -> [e,h,g,f,a,d,c,b]
            assign(y,  vblend_8x32f(u3, u2, [2,1,1,1,2,1,1,1]))))),   # -> [x11,h,g,f,e,d,c,b]

));

#==============================================================================================

SIMD_ISA_DB.addISA(AVX_4x64f);
SIMD_ISA_DB.addISA(AVX_8x32f);


text:
#############################################################################
####
##
#W  anupqprop.gi               ANUPQ package                    Werner Nickel
#W                                                                Greg Gamble
##
##  Installs methods for properties and attributes needed for ANUPQ.
##
#Y  Copyright (C) 2001  Lehrstuhl D fuer Mathematik,  RWTH Aachen,  Germany
##

#############################################################################
##
#F  SET_PQ_PROPS_AND_ATTRS( <G>, <func> )
##
InstallGlobalFunction( SET_PQ_PROPS_AND_ATTRS, function( G, func )
local S;
  if not( HasIsPGroup(G) and IsPGroup(G) ) then
    Error( "supplied group is not known to be a p-group\n" );
  fi;
  S := PqStandardPresentation(G : Prime := PrimePGroup(G));
  SetIsCapable( G, IsCapable(S) );
  SetNuclearRank( G, NuclearRank(S) );
  SetMultiplicatorRank( G, MultiplicatorRank(S) );
  return func(G);
end );

#############################################################################
##
#M  IsCapable( <G> )
##
InstallMethod( IsCapable, "fp p-groups", true, [ IsFpGroup ], 0,
  G -> SET_PQ_PROPS_AND_ATTRS(G, IsCapable)
);

InstallMethod( IsCapable, "pc p-groups", true, [ IsPcGroup ], 0,
  G -> SET_PQ_PROPS_AND_ATTRS(G, IsCapable)
);

#############################################################################
##
#M  NuclearRank( <G> )
##
InstallMethod( NuclearRank, "fp p-groups", [ IsFpGroup ], 0,
  G -> SET_PQ_PROPS_AND_ATTRS(G, NuclearRank)
);

InstallMethod( NuclearRank, "pc p-groups", [ IsPcGroup ], 0,
  G -> SET_PQ_PROPS_AND_ATTRS(G, NuclearRank)
);

#############################################################################
##
#M  MultiplicatorRank( <G> )
##
InstallMethod( MultiplicatorRank, "fp p-groups", [ IsFpGroup ], 0,
  G -> SET_PQ_PROPS_AND_ATTRS(G, MultiplicatorRank)
);

InstallMethod( MultiplicatorRank, "pc p-groups", [ IsPcGroup ], 0,
  G -> SET_PQ_PROPS_AND_ATTRS(G, MultiplicatorRank)
);

#E  anupqprop.gi  . . . . . . . . . . . . . . . . . . . . . . . . . ends here


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

Class(BaseVGath, BaseMat, SumsBase, rec(
    isVGath := true,
    sums := self >> self,
    isReal := self >> true,
    dims := self >> Error(ObjId(self), ".dims() is undefined"),
    needInterleavedLeft := False,
    needInterleavedRight := False,
    conjTranspose := self >> self.transpose(),
));

Class(BaseSVGath, BaseMat, SumsBase, rec(
    isVGath := true,
    isSVGath := true,
    #-----------------------------------------------------------------------
    abbrevs := [
	(func, v, sv) -> Checked(
	    IsFunction(func) or IsFuncExp(func), IsInt(v), IsInt(sv),
	    [func, v, sv, Cond(CanBeFullyEvaluated(func.domain()), EvalScalar((func.domain()*sv mod v) / sv), Unk(TInt))]),
	(func, v, sv, rem) -> Checked(
	    IsFunction(func) or IsFuncExp(func), IsInt(v), IsInt(sv), rem _is Unk or IsInt(_unwrap(rem)),
	    [func, v, sv, rem]),
    ],
    #-----------------------------------------------------------------------
    new := (self, func, v, sv, rem) >> SPL(WithBases(self,
        rec(func := func, v := v, sv := sv, rem := rem))).setDims(),
    #-----------------------------------------------------------------------
    dims := self >> Error(ObjId(self), ".dims() is undefined"),
    #-----------------------------------------------------------------------
    print := (self,i,is) >> Print(self.__name__, "(", self.func, ", ",
        self.v, ", ", self.sv, ", ", self.rem, ")", self.printA()),
    #-----------------------------------------------------------------------
    getConstantRem := self >> Cond(
	self.rem _is Unk and not IsSymbolic(EvalScalar(self.func.domain())),
	    EvalScalar(imod(self.func.domain()*self.sv, self.v) / self.sv),
	self.rem _is Unk,
	    Error("Can't generate sub-vector gather/scatter code for <self>, because its dimensions are unknown"),
	self.rem
    ),
    #-----------------------------------------------------------------------
    conjTranspose := self >> self.transpose(),

    sums := self >> self,
    isReal := self >> true,
    needInterleavedLeft := False,
    needInterleavedRight := True
));

IsVGath := x -> IsRec(x) and IsBound(x.isVGath) and x.isVGath;
IsSVGath := x -> IsRec(x) and IsBound(x.isSVGath) and x.isSVGath;

#F ==========================================================================
#F VGath(<func>, <v>) - vector gather (read) matrix, assumes aligned input
#F
#F <func> must be a vector index function returning vector indices
#F
#F VGath(f, v) == Gath(fTensor(f, fId(v)))
#F ==========================================================================
Class(VGath, BaseVGath, SumsBase, rec(
    #-----------------------------------------------------------------------
    #F NOTE: rChildren: does not expose .v!
    rChildren := self >> [self.func],
    rSetChild := rSetChildFields("func"),
    from_rChildren := (self, rch) >> ObjId(self)(rch[1], self.v),
    #-----------------------------------------------------------------------
    new := (self, func, v) >> SPL(WithBases(self,
        rec(func := func, v := v))).setDims(),
    #-----------------------------------------------------------------------
    dims := self >> [self.v*self.func.domain(), self.v*self.func.range()],
    #-----------------------------------------------------------------------
    area := self >> self.func.domain() * self.v,
    transpose := self >> VScat(self.func, self.v),
    #-----------------------------------------------------------------------
    print := (self,i,is) >> Print(self.name, "(", self.func, ", ", self.v,")", self.printA()),
    #-----------------------------------------------------------------------
    toAMat := self >> let(v:=self.v, n := EvalScalar(v*self.func.domain()),
        N := EvalScalar(v*self.func.range()),
        func := fTensor(self.func, fId(v)).lambda(),
        AMatMat(List([0..n-1], row -> BasisVec(N, EvalScalar(func.at(row).ev()))))),
));

#F ==========================================================================
#F VGath_u(<func>, <v>) -- assumes input is unaligned, and uses scalar <func>
#F
#F VGath_u performs an unaligned gather operations. Thus <func> must be an
#F index function returning scalar, not vector, addresses (unlike VGath).
#F
#F In particular:
#F   VGath(func, vlen) = VGath_u(fTensor(func, fBase(vlen, 0)), vlen)
#F ==========================================================================
Class(VGath_u, VGath, rec(
    new := (self, func, v) >> SPL(WithBases(self,
        rec(func := func, v := v))).setDims(),
    #-----------------------------------------------------------------------
    dims := self >> [self.v*self.func.domain(), self.func.range()],
    #-----------------------------------------------------------------------
    toAMat := self >> let(v:=self.v,
        n := EvalScalar(v*self.func.domain()),
        N := EvalScalar(self.func.range()),
        func := self.func.lambda(),
        AMatMat(List([0..n-1], row -> BasisVec(N, EvalScalar(add(func.at(idiv(row,v)), imod(row, v)).ev()))))),
    #-----------------------------------------------------------------------
    transpose := self >> VScat_u(self.func, self.v)
));

#F ==========================================================================
#F VGath_zero(N, n, <v>) - vector gather matrix
#F ==========================================================================
Class(VGath_zero, BaseVGath, rec(
    #-----------------------------------------------------------------------
    rChildren := self >> [self.N, self.n, self.v],
    rSetChild := rSetChildFields("N", "n", "v"),
    #-----------------------------------------------------------------------
    new := (self, N, n, v) >> SPL(WithBases(self,
        rec(n := n, N := N, v := v))).setDims(),
    #-----------------------------------------------------------------------
    dims := self >> [self.v * self.n, self.v * self.N],
    #-----------------------------------------------------------------------
    transpose := self >> VScat_zero(self.N, self.n, self.v),
    toAMat := self >> VGath(fAdd(self.N, self.n, 0), self.v).toAMat(),
 ));

#F ==========================================================================
#F VGath_dup(<func>, <v>) - gathers scalar values and replicates <v> times to
#F                          obtain a vector values (read) matrix
#F Example:
#F spiral> PrintMat(MatSPL(VGath_dup(fId(4), 2)));
#F [ [ 1,  ,  ,   ],
#F   [ 1,  ,  ,   ],
#F   [  , 1,  ,   ],
#F   [  , 1,  ,   ],
#F   [  ,  , 1,   ],
#F   [  ,  , 1,   ],
#F   [  ,  ,  , 1 ],
#F   [  ,  ,  , 1 ] ]
#F
#F This is performed using shuffles in SSE/SSE2, or movddup in SSE4.
#F
#F The transpose of VGath_dup is VScat_red (red = reduce) incurs
#F arithmetic ops.
#F
#F spiral> PrintMat(MatSPL(VScat_red(fId(4), 2)));
#F [ [ 1, 1,  ,  ,  ,  ,  ,   ],
#F   [  ,  , 1, 1,  ,  ,  ,   ],
#F   [  ,  ,  ,  , 1, 1,  ,   ],
#F   [  ,  ,  ,  ,  ,  , 1, 1 ] ]
#F
#F Some architectures have special instructions for this operation.
#F For example, SSE4 has 'hadd' == horizontal add.
#F
#F ==========================================================================
Declare(VGath_dup);
Class(VGath_dup, VGath, rec(
    #-----------------------------------------------------------------------
    new := (self, func, v) >> SPL(WithBases(self,
        rec(func := func, v := v))).setDims(),
    #-----------------------------------------------------------------------
    dims := self >> [self.v*self.func.domain(), self.func.range()],
    #-----------------------------------------------------------------------
    transpose := self >> Error("not yet supported"), #VScat_red(self.func, self.v),
    #-----------------------------------------------------------------------
    toAMat := self >> Tensor(
        Gath(self.func),
        Mat(TransposedMat([Replicate(self.v, 1)]))).toAMat(),
    #-----------------------------------------------------------------------
    toloop := (self, bksize) >> let(i := Ind(self.func.domain() / bksize),
	ISum(i, i.range,
            VScat(fTensor(fBase(i), fId(bksize)),self.v) *
            VGath_dup(fCompose(self.func, fTensor(fBase(i), fId(bksize))),self.v)))
    #-----------------------------------------------------------------------
));

#F ==========================================================================
#F VGath_dup_x_I(<func>, <v>) - gathers vector values and replicates <v> times
#F NOTE: EXPLAIN THIS
#F ==========================================================================
VGath_dup_x_I := (func, v) -> Cond(
    func.domain() = 1,
        VReplicate(v),
    # else
    let(i := Ind(func.domain()),
	ISum(i, i.range,
	    VScat(fTensor(fBase(i),fId(v)), v) *
	    VReplicate(v) *
	    VGath(fBase(i),v)) *
        VGath(func,v)
    )
);

#F ==========================================================================
#F VGath_sv(<func>, <v>, <sv>) - vector gather (read) matrix on subvectors
#F NOTE: EXPLAIN THIS
#F ==========================================================================
Class(VGath_sv, BaseSVGath, rec(
    #-----------------------------------------------------------------------
    rChildren := self >> [self.func],    # self >> [self.func, self.v, self.sv],
    rSetChild := rSetChildFields("func"), # rSetChildFields("func", "v", "sv"),
    from_rChildren := (self, rch) >> ObjId(self)(rch[1], self.v, self.sv, self.rem),
    #-----------------------------------------------------------------------
    dims := self >> [_roundup(self.sv * self.func.domain(), self.v),
                     self.sv * self.func.range() ],
    #-----------------------------------------------------------------------
    area := self >> 2 * self.func.domain() * self.sv,
    #-----------------------------------------------------------------------
    transpose := self >> VScat_sv(self.func, self.v, self.sv, self.rem),
    #-----------------------------------------------------------------------
    toAMat := self >> let(sv:=self.sv,
        n := EvalScalar(sv * self.func.domain()), N := EvalScalar(Cols(self)),
        nv := EvalScalar(Rows(self)),
	func := fTensor(self.func, fId(sv)).lambda(),
        Scat(fAdd(nv, n, 0)).toAMat() *
        AMatMat(List([0..n-1], row -> BasisVec(N, EvalScalar(func.at(row)))))),
));

#F ==========================================================================
#F RCVGath_sv(<func>, <v>, <sv>) - vector gather (read) matrix on subvectors
#F NOTE: EXPLAIN THIS
#F ==========================================================================
Class(RCVGath_sv, BaseSVGath, rec(
    #-----------------------------------------------------------------------
    rChildren := self >> [self.func],    # self >> [self.func, self.v, self.sv],
    rSetChild := rSetChildFields("func"), # rSetChildFields("func", "v", "sv"),
    from_rChildren := (self, rch) >> ObjId(self)(rch[1], self.v, self.sv, self.rem),
    #-----------------------------------------------------------------------
    dims := self >> let(
        n := self.func.domain() * self.sv, nv := _roundup(n, self.v),
	N := self.func.range() * self.sv,
        2 * [nv, N]),
    #-----------------------------------------------------------------------
    area := self >> 2 * self.func.domain() * self.sv,
    #-----------------------------------------------------------------------
    transpose := self >> RCVScat_sv(self.func, self.v, self.sv, self.rem),
    #-----------------------------------------------------------------------
    toAMat := self >> let(
	sv := EvalScalar(self.sv), n := EvalScalar(self.func.domain()),
	N  := EvalScalar(self.func.range()),
        func := fTensor(self.func, fId(sv), fId(2)).lambda(),
        Scat(fTensor(fAdd(_roundup(n, self.v), n, 0), fId(2))).toAMat() *
	AMatMat(List([0..2*n-1], row -> BasisVec(2*N, EvalScalar(func.at(row).ev()))))),
));

#F ==========================================================================
#F VStretchGath(<func>, <part>, <v>) -
#F   vector gather (read) matrix on subvectors
#F   reads according to func but fills up entire vectors
#F NOTE: EXPLAIN THIS
#F ==========================================================================
Class(VStretchGath, BaseSVGath, rec(
    #-----------------------------------------------------------------------
    rChildren := self >> [self.func],
    rSetChild := rSetChildFields("func"),
    from_rChildren := (self, rch) >> ObjId(self)(rch[1], self.part, self.v),
    #-----------------------------------------------------------------------
    dims := self >> let(
	n  := self.func.domain(), N := self.func.range(),
	nv := self.part * _roundup(n/self.part, self.v),
        [nv, N]),
    #-----------------------------------------------------------------------
    abbrevs := [],
    new := (self, func, part, v) >> SPL(WithBases(self,
        rec(func := func, part := part, v := v))).setDims(),
    #-----------------------------------------------------------------------
    area := self >> self.func.domain(),
    #-----------------------------------------------------------------------
    transpose := self >> VStretchScat(self.func, self.part, self.v),
    #-----------------------------------------------------------------------
    print := (self,i,is) >> Print(self.name, "(", self.func, ", ",
        self.part, ", ", self.v, ")", self.printA()),
    #-----------------------------------------------------------------------
    toAMat := self >> let(
	n  := EvalScalar(self.func.domain()),
	N  := EvalScalar(self.func.range()),
	nv := self.part * _roundup(n/self.part, self.v),
        Tensor(I(self.part), Scat(fAdd(nv/self.part, n/self.part, 0))).toAMat() *
	AMatMat(List([0..n-1], row ->
		    BasisVec(N, EvalScalar(self.func.at(row).ev()))))),
));

#F ==========================================================================
#F vRCStretchGath(<func>, <parts>, <v>) -
#F   vector gather (read) matrix on subvectors
#F   reads according to func but fills up entire vectors
#F NOTE: EXPLAIN THIS
#F ==========================================================================
Class(vRCStretchGath, BaseSVGath, rec(
    #-----------------------------------------------------------------------
    rChildren := self >> [self.func],
    rSetChild := rSetChildFields("func"),
    from_rChildren := (self, rch) >> ObjId(self)(rch[1], self.part, self.v),
    #-----------------------------------------------------------------------
    dims := self >> let(
	n  := self.func.domain(), N := self.func.range(),
	nv := self.part * _roundup(n/self.part, self.v/2),
	2*[nv, N]),
    #-----------------------------------------------------------------------
    abbrevs := [],
    new := (self, func, part, v) >> SPL(WithBases(self,
        rec(func := func, part := part, v := v))).setDims(),
    #-----------------------------------------------------------------------
    area := self >> 2 * self.func.domain(),
    #-----------------------------------------------------------------------
    transpose := self >> vRCStretchScat(self.func, self.part, self.v),
    #-----------------------------------------------------------------------
    print := (self,i,is) >> Print(self.name, "(", self.func, ", ",
        self.part, ", ", self.v, ")", self.printA()),
    #-----------------------------------------------------------------------
    toAMat := self >> let(
	n  := self.func.domain(), N := self.func.range(),
	nv := self.part * _roundup(n/self.part, self.v/2),
        func := fTensor(self.func, fId(2)),
        Tensor(I(self.part),
	       Scat(fAdd(nv/self.part, n/self.part, 0)), I(2)).toAMat() *
	AMatMat(List([0..2*n-1], row -> BasisVec(2*N, EvalScalar(func.at(row).ev()))))),
));

#F ==========================================================================
#F RCVStretchGath(<func>, <parts>, <v>) -
#F   vector gather (read) matrix on subvectors
#F   reads according to func but fills up entire vectors
#F NOTE: EXPLAIN THIS
#F ==========================================================================
Class(RCVStretchGath, BaseSVGath, rec(
    #-----------------------------------------------------------------------
    rChildren := self >> [self.func],
    rSetChild := rSetChildFields("func"),
    from_rChildren := (self, rch) >> ObjId(self)(rch[1], self.part, self.v),
    #-----------------------------------------------------------------------
    dims := self >> let(
	n  := self.func.domain(), N := self.func.range(),
	nv := self.part * _roundup(n/self.part, self.v),
	2*[nv, N]),
    #-----------------------------------------------------------------------
    abbrevs := [],
    new := (self, func, part, v) >> SPL(WithBases(self,
        rec(func := func, part := part, v := v))).setDims(),
    #-----------------------------------------------------------------------
    area := self >> 2 * self.func.domain() * self.v,
    #-----------------------------------------------------------------------
    transpose := self >> RCVStretchScat(self.func, self.part, self.v),
    #-----------------------------------------------------------------------
    print := (self,i,is) >> Print(self.name, "(", self.func, ", ",
        self.part, ", ", self.v, ")", self.printA()),
    #-----------------------------------------------------------------------
    toAMat := self >> Tensor(VStretchGath(self.func, self.part, self.v),
                             I(2)).toAMat(),
));

#F ==========================================================================
#F VGath_pc(N, n, ofs, v) - gather for unaligned partial contiguous vectors
#F   reads <n> contiguous values out of <N> at an offset <ofs>
#F   outputs ceil(<n> / <v>) vectors
#F
#F   .N   -- input size, only affects dimensions of VGath_pc matrix
#F   .n   -- how many contiguous points we are reading
#F   .ofs -- starting offset
#F   .v   -- vector length
#F
#F Example:
#F
#F spiral> PrintMat(MatSPL(VGath_pc(5, 3, 0, 4)));
#F [ [ 1,  ,  ,  ,   ],
#F   [  , 1,  ,  ,   ],
#F   [  ,  , 1,  ,   ],
#F   [  ,  ,  ,  ,   ] ]
#F
#F ==========================================================================
Class(VGath_pc, BaseSVGath, rec(
    #-----------------------------------------------------------------------
    rChildren := self >> [self.N, self.n, self.ofs, self.v],
    rSetChild := rSetChildFields("N", "n", "ofs", "v"),
    #-----------------------------------------------------------------------
    abbrevs := [],
    new := (self, N, n, ofs, v) >> SPL(WithBases(self, rec(
        N := N, n := n, ofs := ofs, v := v))).setDims(),
    #-----------------------------------------------------------------------
    dims := self >> [ _roundup(self.n, self.v), self.N ],
    #-----------------------------------------------------------------------
    area := self >> self.n,
    #-----------------------------------------------------------------------
    transpose := self >> VScat_pc(self.N, self.n, self.ofs, self.v),
    #-----------------------------------------------------------------------
    toAMat := self >> let(n := self.n, N := self.N, nv := _roundup(n, self.v),
        Scat(fAdd(nv, n, 0)).toAMat() * Gath(fAdd(N, n, self.ofs)).toAMat()),
    #-----------------------------------------------------------------------
    print :=  (self, i, is) >> Print(self.__name__, "(", self.N,
        ", ", self.n, ", ", self.ofs, ", ", self.v, ")", self.printA()),
));

#F ==========================================================================
#F IxVGath_pc(k, N, n, ofs, v)
#F
#F Performs multiple stacked VGath_pc operations
#F
#F IxVGath_pc(k, N, n, ofs, v) = I(k) tensor VGath_pc(N, n, ofs, v)
#F
#F See also: VGath_pc
#F
#F Example:
#F
#F spiral> PrintMat(MatSPL(VGath_pc(5, 3, 0, 4)));
#F [ [ 1,  ,  ,  ,   ],
#F   [  , 1,  ,  ,   ],
#F   [  ,  , 1,  ,   ],
#F   [  ,  ,  ,  ,   ] ]
#F
#F spiral> PrintMat(MatSPL(IxVGath_pc(2, 5, 3, 0, 4)));
#F [ [ 1,  ,  ,  ,  ,  ,  ,  ,  ,   ],
#F   [  , 1,  ,  ,  ,  ,  ,  ,  ,   ],
#F   [  ,  , 1,  ,  ,  ,  ,  ,  ,   ],
#F   [  ,  ,  ,  ,  ,  ,  ,  ,  ,   ],
#F   [  ,  ,  ,  ,  , 1,  ,  ,  ,   ],
#F   [  ,  ,  ,  ,  ,  , 1,  ,  ,   ],
#F   [  ,  ,  ,  ,  ,  ,  , 1,  ,   ],
#F   [  ,  ,  ,  ,  ,  ,  ,  ,  ,   ] ]
#F
#F ==========================================================================
Class(IxVGath_pc, BaseSVGath, rec(
    #-----------------------------------------------------------------------
    rChildren := self >> [self.k, self.N, self.n, self.ofs, self.v],
    rSetChild := rSetChildFields("k", "N", "n", "ofs", "v"),
    #-----------------------------------------------------------------------
    abbrevs := [],
    new := (self, k, N, n, ofs, v) >> SPL(WithBases(self, rec(
        k := k, N := N, n:= n, ofs := ofs, v := v))).setDims(),
    #-----------------------------------------------------------------------
    dims := self >> [ self.k*_roundup(self.n, self.v), self.k*self.N ],
    #-----------------------------------------------------------------------
    area := self >> self.k * self.n,
    #-----------------------------------------------------------------------
    transpose := self >> IxVScat_pc(self.k, self.N, self.n, self.ofs, self.v),
    #-----------------------------------------------------------------------
    toloop := (self, bksize) >> let(
	j := Ind(self.k),
	gath := VGath_pc(self.k*self.N, self.n, self.ofs + j*self.N, self.v),
	nvecs := Rows(gath) / self.v,
	ISum(j, j.range,
             VScat(HH(self.k*nvecs, nvecs, j*nvecs, [1]), self.v) *
	     gath).split(bksize)),
    #-----------------------------------------------------------------------
    toAMat := self >>
        Tensor(I(self.k), VGath_pc(self.N, self.n, self.ofs, self.v)).toAMat(),
    #-----------------------------------------------------------------------
    print :=  (self, i, is) >> Print(self.__name__, "(", self.k, ", ", self.N,
        ", ", self.n, ", ", self.ofs, ", ", self.v, ")", self.printA()),
));

#F ==========================================================================
#F IxRCVGath_pc(k, N, n, ofs, v)
#F   vector gather (read) matrix reading unaligned
#F   partial contiguous COMPLEX vectors
#F
#F   IxRCVGath_pc(..params..) == IxVGath_pc(..params..) tensor I(2)
#F
#F  Note: IxRCVGath_pc(k, N, n, ofs, v) == IxVGath_pc(k, 2*N, 2*n, 2*ofs, 2*v)
#F        observe 2*v inside IxVGath_pc.
#F
#F        IxVGath_pc(j, 2*N, 2*n, 2*ofs, v), however, is not exactly equal
#F        it has a different # of rows, if n < v/2
#F ==========================================================================
Class(IxRCVGath_pc, BaseSVGath, rec(
    #-----------------------------------------------------------------------
    rChildren := self >> [self.k, self.N, self.n, self.ofs, self.v],
    rSetChild := rSetChildFields("k", "N", "n", "ofs", "v"),
    #-----------------------------------------------------------------------
    abbrevs := [],
    new := (self, k, N, n, ofs, v) >> SPL(WithBases(self, rec(
        k := k, N := N, n := n, ofs := ofs, v := v))).setDims(),
    #-----------------------------------------------------------------------
    dims := self >> [ 2*self.k*_roundup(self.n, self.v),  2*self.k*self.N ],
    #-----------------------------------------------------------------------
    area := self >> 2 * self.k * self.n,
    #-----------------------------------------------------------------------
    transpose := self >> IxRCVScat_pc(self.k, self.N, self.n, self.ofs, self.v),
    #-----------------------------------------------------------------------
    toAMat := self >>
        Tensor(IxVGath_pc(self.k, self.N, self.n, self.ofs, self.v), I(2)).toAMat(),
    #-----------------------------------------------------------------------
    print :=  (self, i, is) >> Print(self.__name__, "(", self.k, ", ", self.N,
        ", ", self.n, ", ", self.ofs, ", ", self.v, ")", self.printA()),
));


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

_computeExpType := function(e)
    local t;
    if IsBound(e.computeType) then t := e.computeType();
#    elif Length(e.args)=0 then t := TUnknown;
    else t := UnifyTypes(List(e.args, x->x.t));
    fi;
    e.t := t;
    return e;
end;

PropagateTypes := c -> SubstBottomUp(c,
    @.cond(e->IsExp(e) and not IsBound(e.t)),
    _computeExpType);


text:
InstallGlobalFunction(MitM_GAPToXML,
obj -> MitM_OMRecToXML(MitM_GAPToOMRec(obj)));

InstallGlobalFunction(MitM_XMLToGAP,
str -> MitM_OMRecToGAP(MitM_XMLToOMRec(str)));

InstallGlobalFunction(MitM_RoundTripGAP,
obj -> MitM_XMLToGAP(MitM_GAPToXML(obj)));

InstallGlobalFunction(MitM_RoundTripXML,
function(str)
    local r;
    r := MitM_XMLToGAP(str);
    if r.success <> true then
        return r;
    fi;
    return MitM_GAPToXML(r.result);
end);

InstallGlobalFunction(MitM_Print,
function(obj)
    Print(MitM_OMRecToXML(MitM_OMRecToOMOBJRec(MitM_GAPToOMRec(obj))), "\n");
end);

InstallGlobalFunction(MitM_OMRecToOMOBJRec,
r -> OMOBJ(r));

max_stars_repo_name: markuspf/MathsInTheMiddle

max_issues_repo_name: markuspf/MathsInTheMiddle

max_forks_repo_name: markuspf/MathsInTheMiddle


text:
IsExistingFile("input.txt");
IsDirectoryPath("docs");
IsExistingFile("/input.txt");
IsDirectoryPath("/docs");

ext: gap

max_stars_repo_name: mullikine/RosettaCodeData

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData