text:
--
-- An LPG Lexer Template Using lpg.jar
--
-- An instance of this template must have a %Export section and the export_terminals option
-- There must be only one non-terminal, the start symbol, for the keywords
-- The action for each keyword should be a call to %setResult(terminal_symbol)
--
-- Macro that may be redefined in an instance of this template
--
--     %eof_char
--
-- B E G I N N I N G   O F   T E M P L A T E   KeywordTemplateF (Similar to KeywordTemplateD)
--
%Options programming_Language=typescript,margin=4
%Options table
%options action-block=("*.ts", "/.", "./")
%options ParseTable=ParseTable
%Options prefix=Char_
%Options single-productions

%Globals
    /.
    import { %prs_type } from ".\/%prs_type";
    import { %sym_type } from ".\/%sym_type";
    import { %exp_type } from ".\/%exp_type";
    ./
%End

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

%Define
    --
    -- Macro that may be respecified in an instance of this template
    --
    $eof_char /.%sym_type%.%prefix%EOF%suffix%./

    --
    -- Macros useful for specifying actions
    --
    $setResult /.this.keywordKind[%rule_number] = ./

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

    $BeginAction /.%Header./

    $EndAction /../

    $BeginJava /.%BeginAction./

    $EndJava /.%EndAction./
%End

%Headers
    /.
    export class %action_type extends %prs_type
    {
        private inputChars : string;
        private   keywordKind  : number[] = new Array(%num_rules + 1);

        public  getKeywordKinds() : number[] { return this.keywordKind; }

        public  lexer(curtok : number, lasttok : number) : number
        {
            let current_kind = %action_type.getKind(this.inputChars.charCodeAt(curtok)),
                act;

            for (act = this.tAction(this.START_STATE, current_kind);
                 act > this.NUM_RULES && act < this.ACCEPT_ACTION;
                 act = this.tAction(act, current_kind))
            {
                curtok++;
                current_kind = (curtok > lasttok
                                       ? %eof_char
                                       : %action_type.getKind(this.inputChars.charCodeAt(curtok)));
            }

            if (act > this.ERROR_ACTION)
            {
                curtok++;
                act -= this.ERROR_ACTION;
            }

            return this.keywordKind[act == this.ERROR_ACTION  || curtok <= lasttok ? 0 : act];
        }

        public setInputChars(inputChars : string ) : void  { this.inputChars = inputChars; }

    ./
%End

%Rules
    /.

        constructor( inputChars : string,  identifierKind : number)
        {
            super();
            this.inputChars = inputChars;
            this.keywordKind[0] = identifierKind;
    ./
%End

%Trailers
    /.
            for (let i : number = 0; i < this.keywordKind.length; i++)
            {
                if (this.keywordKind[i] == 0)
                    this.keywordKind[i] = identifierKind;
            }
        }
    }
    ./
%End

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


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

###############################################################################
##
#R  IsTreeHomomorphismRep
##
DeclareRepresentation("IsTreeHomomorphismRep",
                      IsComponentObjectRep and IsAttributeStoringRep,
                      ["states", "perm", "deg"]);

###############################################################################
##
#R  IsTreeHomomorphismFamilyRep
##
DeclareRepresentation("IsTreeHomomorphismFamilyRep",
                      IsComponentObjectRep and IsAttributeStoringRep,
                      ["spher_index", "top_deg"]);

###############################################################################
##
##  AG_CreatedTreeHomomorphismFamilies
##
##  Contains all created TreeHomomorphismFamily objects; for each spherical
##  index there exists one family, to which all objects created with TreeHomomorphism
##  belong.
##
BindGlobal("AG_CreatedTreeHomomorphismFamilies", rec(ind := [], fam := []));

###############################################################################
##
#M  TreeHomomorphismFamily(<sph_ind>)
##
InstallMethod(TreeHomomorphismFamily, [IsRecord],
function(sph_ind)
  local fam, pos;

  sph_ind := AG_ReducedSphericalIndex(sph_ind);
  if sph_ind in AG_CreatedTreeHomomorphismFamilies.ind then
    for fam in AG_CreatedTreeHomomorphismFamilies.fam do
      if fam!.spher_index = sph_ind then
        return fam;
      fi;
    od;
  fi;

  fam := NewFamily(Concatenation("Automorphisms of ", sph_ind.start, ", (", sph_ind.period, ")-tree"),
                   IsTreeHomomorphism, IsTreeHomomorphism,
                   IsTreeHomomorphismFamily and IsTreeHomomorphismFamilyRep);
  fam!.spher_index := sph_ind;
  fam!.top_deg := AG_TopDegreeInSphericalIndex(sph_ind);

  AddSet(AG_CreatedTreeHomomorphismFamilies.ind, sph_ind);
  Add(AG_CreatedTreeHomomorphismFamilies.fam, fam);

  return fam;
end);

###############################################################################
##
#M  TreeHomomorphism (<states>, <tr>)
##
InstallMethod(TreeHomomorphism,
              [IsList and IsTreeHomomorphismCollection, IsObject],
function(states, perm)
  local top_deg, bot_deg, ind, fam, a;

  if not IsPerm(perm) and not IsTransformation(perm) then
    Error("The second argument ",perm, "must be a permutation or transformation");
  fi;

  if perm^-1<>fail and
     ForAll(states, IsTreeAutomorphism)
  then
    return TreeAutomorphism(states, AG_PermFromTransformation(perm));
  fi;

  top_deg := Length(states);

  if IsPerm(perm) then
    if not IsOne(perm) and top_deg < Maximum(MovedPoints(perm)) then
      Error("The root permutation ", perm, " must move only points from 1 to the degree ", top_deg, " of the tree");
    fi;
  else
    if not IsOne(perm) and top_deg < DegreeOfTransformation(perm) then
      Error("The root transformation ", perm, " must move only points from 1 to the degree ", top_deg, " of the tree");
    fi;
  fi;

  bot_deg := DegreeOfTree(states[1]);
  ind := rec(start := [top_deg], period := [bot_deg]);
  fam := TreeHomomorphismFamily(ind);

  return Objectify(NewType(fam, IsTreeHomomorphism and IsTreeHomomorphismRep),
                   rec(states := ShallowCopy(states),
                       perm := perm,
                       deg := top_deg));
end);

###############################################################################
##
#M  TreeHomomorphism (<states_list>, <perm>)
##
InstallMethod(TreeHomomorphism, [IsList, IsTransformation],
function(states, perm)
  local autom, nstates, s;

  autom := fail;

  for s in states do
    if IsTreeHomomorphism(s) then
      autom := s;
      break;
    elif not IsOne(s) then
      Error("Invalid state `", s, "'");
    fi;
  od;

  if autom = fail then
    Error("Can't create an automaton with all trivial states ",
          "without information about the tree");
  fi;

  nstates := List(states, function(s)
                            if IsOne(s) then
                              return One(autom);
                            else
                              return s;
                            fi;
                          end);

  return TreeHomomorphism(nstates, perm);
end);

###############################################################################
##
#M  TreeHomomorphism(<state_1>, <state_2>, ..., <state_n>, <perm>)
##
InstallMethod(TreeHomomorphism, [IsObject, IsObject, IsTransformation],
  function(a1, a2, perm) return TreeHomomorphism([a1, a2], perm); end);
InstallMethod(TreeHomomorphism, [IsObject, IsObject, IsObject, IsTransformation],
  function(a1, a2, a3, perm) return TreeHomomorphism([a1, a2, a3], perm); end);
InstallMethod(TreeHomomorphism, [IsObject, IsObject, IsObject, IsObject, IsTransformation],
  function(a1, a2, a3, a4, perm) return TreeHomomorphism([a1, a2, a3, a4], perm); end);
InstallMethod(TreeHomomorphism, [IsObject, IsObject, IsPerm],
  function(a1, a2, perm) return TreeHomomorphism([a1, a2], perm); end);
InstallMethod(TreeHomomorphism, [IsObject, IsObject, IsObject, IsPerm],
  function(a1, a2, a3, perm) return TreeHomomorphism([a1, a2, a3], perm); end);
InstallMethod(TreeHomomorphism, [IsObject, IsObject, IsObject, IsObject, IsPerm],
  function(a1, a2, a3, a4, perm) return TreeHomomorphism([a1, a2, a3, a4], perm); end);

###############################################################################
##
#M  ViewObj(<a>)
##
InstallMethod(ViewObj, [IsTreeHomomorphism],
function (a)
    local deg, printword, i, perm, states;

    states := Sections(a);
    deg := Length(states);
    perm := TransformationOnLevel(a, 1);

    Print("(");
    for i in [1..deg] do
        View(states[i]);
        if i <> deg then Print(", "); fi;
    od;
    Print(")");
    if not IsOne(perm) then
      AG_PrintTransformation(perm);
    fi;
end);

###############################################################################
##
#M  PrintObj(<a>)
##
InstallMethod(PrintObj, "for [IsTreeHomomorphism and IsTreeHomomorphismRep]",
                             [IsTreeHomomorphism and IsTreeHomomorphismRep],
function (a)
    local deg, i, states, perm;

    states := Sections(a);
    deg := Length(states);
    perm := TransformationOnLevel(a, 1);

    Print("(");
    for i in [1..deg] do
      if IsAutom(a!.states[i]) then
        View(a!.states[i]);
      else
        Print(a!.states[i]);
      fi;
      if i <> deg then Print(", "); fi;
    od;
    Print(")");
    if not IsOne(perm) then
      AG_PrintTransformation(perm);
    fi;
end);

###############################################################################
##
#M  String(<a>)
##
InstallMethod(String, "for [IsTreeHomomorphism]", [IsTreeHomomorphism],
function (a)
    local deg, printword, i, perm, states, str;

    states := Sections(a);
    deg := Length(states);
    perm := TransformationOnLevel(a, 1);
    str:= "(";

    for i in [1..deg] do
        Append(str, String(states[i]));
        if i <> deg then Append(str, ", "); fi;
    od;
    Append(str, ")");
    if not IsOne(perm) then
      Append(str, AG_TransformationString(perm));
    fi;
    return str;
end);

###############################################################################
##
#M  SphericalIndex (<a>)
##
InstallMethod(SphericalIndex, [IsTreeHomomorphism and IsTreeHomomorphismRep],
function(a)
  return FamilyObj(a)!.spher_index;
end);
InstallMethod(TopDegreeOfTree, [IsTreeHomomorphism and IsTreeHomomorphismRep],
function(a)
  return FamilyObj(a)!.top_deg;
end);

###############################################################################
##
#M  TransformationOnLevel (<a>, <k>)
##
InstallMethod(TransformationOnLevelOp, "for [IsTreeHomomorphism, IsPosInt]",
              [IsTreeHomomorphism, IsPosInt],
function(a, k)
  local states, top, first_level, i, j, d1, d2, permuted, p;

  if k = 1 then
    return TransformationOnFirstLevel(a);
  fi;

  # TODO: it is unnesessarily greedy, it could check whether there
  # are trivial permutations below
  d1 := DegreeOfTree(a);
  d2 := 1;
  for i in [2 .. k] do
    d2 := d2 * DegreeOfLevel(a, i);
  od;
  states := Sections(a);
  top := TransformationOnFirstLevel(a);
  first_level := List(states, s -> TransformationOnLevel(s, k-1));
  permuted := [];
  for i in [1..d1] do
    for j in [1..d2] do
      permuted[d2*(i-1) + j] := d2*(i^top - 1) + j^first_level[i];
    od;
  od;

#  p := PermList(permuted);
#  if p = fail then
#    p := Transformation(permuted);
#  fi;
#  return p;

  return Transformation(permuted);
end);

InstallMethod(TransformationOnFirstLevel, [IsTreeHomomorphism and IsTreeHomomorphismRep],
function(a)
  return AsTransformation(a!.perm);
end);

###############################################################################
##
#M  k ^ a
##
InstallMethod(\^, "for [IsPosInt, IsTreeHomomorphism]", [IsPosInt, IsTreeHomomorphism],
function(k, a)
    return k ^ TransformationOnLevel(a, 1);
end);

###############################################################################
##
#M  seq ^ a
##
InstallMethod(\^, "for [IsList, IsTreeHomomorphism]", [IsList, IsTreeHomomorphism],
function(seq, a)
  if Length(seq) = 0 then
    return [];
  elif Length(seq) = 1 then
    return [seq[1]^TransformationOnLevel(a, 1)];
  else
    return Concatenation([seq[1]^TransformationOnLevel(a, 1)],
                         seq{[2..Length(seq)]}^Section(a, seq[1]));
  fi;
end);

# ###############################################################################
# ##
# #M  FixesLevel(<a>, <k>)
# ##
# InstallMethod(FixesLevel, "for [IsTreeHomomorphism, IsPosInt]",
#               [IsTreeHomomorphism, IsPosInt],
# function(a, k)
#   if HasIsSphericallyTransitive(a) then
#     if IsSphericallyTransitive(a) then
#       return false; fi; fi;
#
#   if IsOne(PermOnLevel(a, k)) then
#     Info(InfoAutomGrp, 3, "IsSphericallyTransitive(a): false");
#     Info(InfoAutomGrp, 3, "  a is not transitive on level", k);
#     Info(InfoAutomGrp, 3, "  a = ", a);
#     SetIsSphericallyTransitive(a, false);
#     return true;
#   else
#     return false;
#   fi;
# end);
#
#
# ###############################################################################
# ##
# #M  FixesVertex(<a>, <v>)
# ##
# InstallOtherMethod(FixesVertex,  "for [IsTreeHomomorphism, IsObject]",
#                    [IsTreeHomomorphism, IsObject],
# function(a, v)
#   if HasIsSphericallyTransitive(a) then
#     if IsSphericallyTransitive(a) then
#       Info(InfoAutomGrp, 3, "FixesVertex(a, v): false");
#       Info(InfoAutomGrp, 3, "  IsSphericallyTransitive(a)");
#       Info(InfoAutomGrp, 3, "  a = ", a);
#       return false;
#     fi;
#   fi;
#
#   if v^a = v then
#     Info(InfoAutomGrp, 3, "IsSphericallyTransitive(a): false");
#     Info(InfoAutomGrp, 3, "  a fixes vertex ", v);
#     Info(InfoAutomGrp, 3, "  a = ", a);
#     SetIsSphericallyTransitive(a, false);
#     return true;
#   else
#     return false;
#   fi;
# end);

###############################################################################
##
#M  Section(<a>, <k>)
##
InstallMethod(Section, [IsTreeHomomorphism, IsPosInt],
function(a, k)
  return Sections(a)[k];
end);

InstallMethod(Section, [IsTreeHomomorphism and IsTreeHomomorphismRep, IsPosInt],
function(a, k)
  return a!.states[k];
end);

###############################################################################
##
#M  Section(<a>, <v>)
##
InstallMethod(Section, [IsTreeHomomorphism, IsList],
function(a, v)
  if Length(v) = 1 then
    return Section(a, v[1]);
  else
    return Section(Section(a, v[1]), v{[2..Length(v)]});
  fi;
end);

###############################################################################
##
#M  Sections(<a>)
##
InstallMethod(Sections, [IsTreeHomomorphism and IsTreeHomomorphismRep],
function(a)
  return a!.states;
end);

###############################################################################
##
#M  Sections(a, k)
##
InstallMethod(Sections, "for [IsTreeHomomorphism, IsPosInt]",
                   [IsTreeHomomorphism, IsPosInt],
function(a, level)
  if level = 1 then
    return Sections(a);
  else
    return Concatenation(List(Sections(a), s -> Sections(s, level-1)));
  fi;
end);

InstallMethod(Sections, "for [IsTreeHomomorphism, IsInt and IsZero]", [IsTreeHomomorphism, IsInt and IsZero],
function(a, level)
  return [a];
end);

###############################################################################
##
#M  Decompose(<a>, <k>)
##
InstallMethod(Decompose, "for [IsTreeHomomorphism, IsPosInt]",
              [IsTreeHomomorphism, IsPosInt],
function(a, level)
  return TreeHomomorphism(Sections(a, level), TransformationOnLevel(a, level));
end);

InstallMethod(Decompose, [IsTreeHomomorphism, IsInt and IsZero],
function(a, level)
  return a;
end);

###############################################################################
##
#M  Decompose(<a>)
##
InstallMethod(Decompose, "for [IsTreeHomomorphism]",
              [IsTreeHomomorphism],
function(a)
  return Decompose(a, 1);
end);

###############################################################################
##
#M  IsOne(<a>)
##
InstallMethod(IsOne, "for [IsTreeHomomorphism]",
              [IsTreeHomomorphism],
function(a)
  local s;

  if not IsOne(TransformationOnLevel(a, 1)) then
    return false;
  fi;

  for s in Sections(a) do
    if not IsOne(s) then
      return false;
    fi;
  od;

  return true;
end);

###############################################################################
##
#M  \=(<a1>, <a2>)
##
# TODO: can lead to infinite recursion
InstallMethod(\=, "for [IsTreeHomomorphism, IsTreeHomomorphism]", ReturnTrue,
              [IsTreeHomomorphism, IsTreeHomomorphism],
function(a1, a2)
  return TransformationOnLevel(a1, 1) = TransformationOnLevel(a2, 1) and
          Sections(a1) = Sections(a2);
end);

###############################################################################
##
#M  \<(<a1>, <a2>)
##
InstallMethod(\<, [IsTreeHomomorphism and IsTreeHomomorphismRep,
                   IsTreeHomomorphism and IsTreeHomomorphismRep],
function(a1, a2)
  return AG_TreeHomomorphismCmp(a1, a2) < 0;
end);

###############################################################################
##
##  AG_TreeHomomorphismCmp(a1, a2)
##
##  Global function to be used from IsTreeAutomomorphism too
##
InstallGlobalFunction(AG_TreeHomomorphismCmp,
function(a1, a2)
  local i, cmp;

  cmp := AG_TrCmp(a1!.perm, a2!.perm, a1!.deg);

  if cmp < 0 then
    return -1;
  elif cmp > 0 then
    return 1;
  fi;

  for i in [1..a1!.deg] do
    if a1!.states[i] < a2!.states[i] then
      return -1;
    elif a1!.states[i] > a2!.states[i] then
      return 1;
    fi;
  od;

  return 0;
end);

###############################################################################
##
#M  OneOp(<a>)
##
InstallMethod(OneOp, [IsTreeHomomorphism and IsTreeHomomorphismRep],
function(a)
  return Objectify(NewType(FamilyObj(a), IsTreeHomomorphism and IsTreeHomomorphismRep),
                   rec(states := List([1..a!.deg], i -> One(a!.states[1])),
                       perm := Transformation([1..a!.deg]),
                       deg := a!.deg));
end);

###############################################################################
##
#M  \*(<a1>, <a2>)
##
InstallMethod(\*, [IsTreeHomomorphism and IsTreeHomomorphismRep,
                   IsTreeHomomorphism and IsTreeHomomorphismRep],
function(a1, a2)
  local a;
  a := Objectify(NewType(FamilyObj(a1), IsTreeHomomorphism and IsTreeHomomorphismRep),
        rec(states := List([1..a1!.deg], i -> a1!.states[i] * a2!.states[i^(a1!.perm)]),
            perm := a1!.perm * a2!.perm,
            deg := a1!.deg));
  SetIsActingOnBinaryTree(a, IsActingOnBinaryTree(a1));
  SetIsActingOnRegularTree(a, IsActingOnRegularTree(a1));
  return a;
end);

###############################################################################
##
#M  \[\](<a1>, <a2>)
##
InstallOtherMethod(\[\], [IsTreeHomomorphism, IsPosInt],
function(a, k)
  return Section(a, k);
end);

###############################################################################
##
#M  Representative( <word>, <fam> )
##
InstallMethod(Representative, "for [IsAssocWord, IsTreeHomomorphismFamily]",
              [IsAssocWord, IsTreeHomomorphismFamily],
function( word, fam )
  if IsAutomFamily( fam ) then return Autom( word, fam );
  elif IsSelfSimFamily( fam ) then return SelfSim( word, fam );
  else Error("the family <fam> must be either IsAutomFamily or IsSelfSimFamily");
  fi;
end);

###############################################################################
##
#M  Representative( <word>, <a> )
##
InstallMethod(Representative, "for [IsAssocWord, IsTreeHomomorphism]",
              [IsAssocWord, IsTreeHomomorphism],
function( word, a )
  local fam;
  fam := FamilyObj(a);
  if IsAutomFamily( fam ) then return Autom( word, fam );
  elif IsSelfSimFamily( fam ) then return SelfSim( word, fam );
  else Error("the homomorphism <a> must be either from IsAutomFamily or from IsSelfSimFamily");
  fi;
end);

# ###############################################################################
# ##
# #M  InverseOp(<a>)
# ##
# InstallMethod(InverseOp, "for [IsTreeHomomorphism and IsTreeHomomorphismRep]",
#               [IsTreeHomomorphism and IsTreeHomomorphismRep],
# function(a)
#   local inv;
#   inv := Objectify(NewType(FamilyObj(a), IsTreeHomomorphism and IsTreeHomomorphismRep),
#             rec(states := List([1..a!.deg], i -> a!.states[i^(a!.perm^-1)]^-1),
#                 perm := a!.perm ^ -1,
#                 deg := a!.deg) );
#   SetIsActingOnBinaryTree(inv, IsActingOnBinaryTree(a));
#   SetIsActingOnRegularTree(inv, IsActingOnRegularTree(a));
#   return inv;
# end);
#
# InstallMethod(InverseOp, "for [IsTreeHomomorphism]", [IsTreeHomomorphism],
# function(a)
#   local states, inv_states, perm;
#   states := Sections(a);
#   perm := Inverse(Perm(a));
#   inv_states := List([1..Length(states)], i -> Inverse(states[i^perm]));
#   return TreeHomomorphism(inv_states, perm);
# end);

#E

ext: gi

max_stars_repo_name: gap-packages/automgrp

max_issues_repo_name: gap-packages/automgrp

max_forks_repo_name: gap-packages/automgrp


text:
ITER_POLY_WARN:=false;
Read(Filename(home_dir,"lib/io.gi"));
Read(Filename(home_dir,"handle.gi"));

F:=GF(3^2);
F_gen:=GeneratorsOfField(F)[1];
aaa:=F_gen^2;

sys:=ChevalleyAdj("F",4,F);
alg:=AlgebraicU(sys);
orbs:=UnipotentClasses(alg,"");
infos:=List(AllClasses(orbs),o->handleClassShort(o));
labels:=List(AllClasses(orbs),x->Label(x));

APR:=ring(alg);
avars:=IndeterminatesOfPolynomialRing(APR);
vals:=ShallowCopy(avars);
pr:=positiveRoots(alg);
generic:=Unipotent(alg,List([1..Length(pr)],i->[i,avars[i]]));

Read(Filename(components_dir,"F4char3.gi"));


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

# ==========================================================================
# Sparse
# ==========================================================================
Class(Sparse, BaseMat, rec(
    new := meth(self, L)
        Constraint(IsList(L));
	if not ForAll(L, t -> IsList(t) and Length(t)=3 and
		              IsInt(t[1]) and IsInt(t[2])) then
	    Error("<L> must be a list of triples [i, j, a_(i,j)]");
	fi;
        return SPL(WithBases( self, rec( element := L,
				       dimensions := [ Maximum(List(L, t->t[1])),
				                       Maximum(List(L, t->t[2])) ]
                       )));
    end,
    #-----------------------------------------------------------------------
    dims := self >> [ Maximum(List(self.element, t->t[1])),
	              Maximum(List(self.element, t->t[2])) ],
    #-----------------------------------------------------------------------
    isPermutation := False,
    #-----------------------------------------------------------------------
    isReal := self >> ForAll(self.element, t -> IsRealNumber(t[3])),
    #-----------------------------------------------------------------------
    toAMat := self >> AMatMat(List(MatSparseSPL(self), r -> List(r, EvalScalar))),
    #-----------------------------------------------------------------------
    transpose := meth(self)  # we use CopyFields to copy all fields of self
        local L, t;
	L := [ ];
	for t in self.element do
	    Add(L, [t[2], t[1], t[3]]);
	od;
	return CopyFields(self, rec(element := L,
		                 dimensions := Reversed(self.dims())));
    end,
    conjTranspose := meth(self)  # we use CopyFields to copy all fields of self
        local L, t;
	L := [ ];
	for t in self.element do
	    Add(L, [t[2], t[1], Global.Conjugate(t[3])]);
	od;
	return CopyFields(self, rec(element := L,
		                 dimensions := Reversed(self.dims())));
    end,
    #-----------------------------------------------------------------------
    arithmeticCost := meth(self, costMul, costAddMul)
        local cost, row, elms;
	cost := costMul(0) - costMul(0); # will work even when costMul(0) <> 0
	for row in [1..self.dimensions[1]] do
	    elms := Filtered(self.element, e -> e[1]=row);
	    if Length(elms) > 0 then
		cost := cost + costMul(elms[1][3])
		             + Sum(elms{[2..Length(elms)]}, e -> costAddMul(e[3]));
	    fi;
	od;
	return cost;
    end
));


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

# First, for Pease, need parameters:
# 1. n
# 2. direction (forward or inverse)
# 3. data width
# 4. twiddle width
# 5. bit reverse in/out
# 6. scaling?
# 7. parallelism (streaming width)
# 8. BRAM budget

# bitrev = 0 --> n/n
# bitrev = 1 --> n/r
# bitrev = 2 --> r/n

# direction = 0 --> DFT(n, -1) (forward)
# direction = 1 --> DFT(n, 1) (inverse)

# test:
# webGen(0, 16, 2, 0, 0, 16, 16, 0, 0, 4, -1, 0, "test");
# output in /tmp/spiral/[PID]/test.v

_splhdl_path := "";  #set in .gaprc

webGen := function(arch, size, radix, direction, d_type, d_width, t_width, bitrev, scaling, width, brams, ip, name)
   local opts, t, fwd, scale_str, bram_str;

   # With radix limited to 16, this will avoid more expensive 2/8 and 8/2 CT splits.
#   DFT_CT.children := (nt) -> Map2(Filtered(DivisorPairs(nt.params[1]), divpair -> (divpair[1] <= 4 and divpair[2] <= 4)), (m,n) -> [ DFT(m, nt.params[2] mod m), DFT(n, nt.params[2] mod n) ]);

 if (ip = 0) then
       paradigms.stream.avoidPermPatent := 1;
   fi;

   if (arch = 0) then
       opts := InitStreamHw();
   else
       opts := InitStreamUnrollHw();
   fi;

   fwd := Cond(direction = 0, -1, 1);
   scale_str := Cond(scaling = 0, " ", " -s ");
   bram_str := Cond(brams = -1, " ", ConcatenationString(" -br ", String(brams)));

    if (_splhdl_path = "") then
	Error("Error: path to SPLHDL compiler not set: paradigms.stream._splhdl_path is not bound.");
    fi;

   opts.profile.makeopts.SPLHDL := _splhdl_path; #"/home/pam/web_backends/splhdl/src/splhdl";
   opts.profile.makeopts.OUTNAME := ConcatenationString(String(name), ".v");
   opts.profile.makeopts.GAP := "gap.spl";
   opts.profile.makeopts.WRAP := ConcatenationString(scale_str, bram_str, " -lic -web");

   if (d_type = 0) then
       opts.profile.makeopts.DATATYPE := ConcatenationString("fix ", StringInt(d_width));
   else
       opts.profile.makeopts.DATATYPE := ConcatenationString("flt ", StringInt(d_width));
   fi;

#       opts.profile.makeopts.TWIDTYPE := ConcatenationString("-fixtw", " ", StringInt(t_width));
#   Print("name: ", name, "\n");

   DFTDR_tSPL_Pease.unroll_its := 1;
   DRDFT_tSPL_Pease.unroll_its := 1;

   if ((size = 4) and (width=4) and (radix=2)) then
       radix := 4;
   fi;

  if (bitrev = 0) then
       t := streamDFT(TRC(DFT(size, fwd)).withTags([AStream(2*width)]), radix, opts);
   else
       if (bitrev = 1) then
           t := streamDFT(TRC(DRDFT(size, fwd, radix)).withTags([AStream(2*width)]), radix, opts);
       else
           t := streamDFT(TRC(DFTDR(size, fwd, radix)).withTags([AStream(2*width)]), radix, opts);
       fi;
   fi;

#   Print("t = ", t, "\n\n\n");

   return CMeasure(t, opts);

end;

webGenDCT := function(size, d_type, d_width, brams, width, ip, name)
   local opts, bram_str, t, width;

   if (ip = 0) then
       paradigms.stream.avoidPermPatent := 1;
   fi;

   opts := InitStreamUnrollHw();
   bram_str := Cond(brams = -1, " ", ConcatenationString(" -br ", String(brams)));

    if (_splhdl_path = "") then
	Error("Error: path to SPLHDL compiler not set: paradigms.stream._splhdl_path is not bound.");
    fi;

   opts.profile.makeopts.SPLHDL := _splhdl_path;
   opts.profile.makeopts.OUTNAME := ConcatenationString(String(name), ".v");
   opts.profile.makeopts.GAP := "gap.spl";
   opts.profile.makeopts.WRAP := ConcatenationString(bram_str, " -lic -webdct");

   if (d_type = 0) then
       opts.profile.makeopts.DATATYPE := ConcatenationString("fix ", StringInt(d_width));
   else
       opts.profile.makeopts.DATATYPE := ConcatenationString("flt ", StringInt(d_width));
   fi;

   t := streamDCTUnroll(size, width);
   return CMeasure(t, opts);
end;


text:
##############################################################
##
##           Helps compute the decomposition of
##           the ramification module and the equivariant
##           degree of a multiple of a simple orbit for
##           the Hurwitz curves X with automorphism group
##           G = PSL(2,q), q a prime.
##
##           Using these computations and Borne's formula, the
##           G-module structure of the RR spaces of equivariant
##           divisors can be determined explicitly.
##
## Reference: David Joyner, Amy Ksir, Roger Vogeler,
##            "Group representations on Riemann-Roch spaces
##            of some Hurwitz curves," preprint, 2006.
##
## 6-9-2006, wdj
##############################################################

## To read this in, use for example:
## Read("/home/wdj/gapfiles/curves/hurwitz_crv_rr_sp.gap");
## To log your results, use for example:
## LogTo("/home/wdj/gapfiles/curves/hurwitz_crv_rr_sp1.log");

hurwitz_primes:=function(n)
# Returns the primes < n s.t. PSL(2,p) is a Hurwitz gp
# Assumes n<1000.
local p,L;
L:=[];
for p in Primes do
 if Int((p+1)/7)=((p+1)/7) and p<n then
   L:=Concatenation(L,[p]);
 fi;
 if Int((p-1)/7)=((p-1)/7) and p<n then
   L:=Concatenation(L,[p]);
 fi;
od;
return L;
end;

hurwitz_prime_powers:=function(n)
## Returns the prime powers q=p^3 < n s.t. PSL(2,q)
## is a Hurwitz gp
## Assumes n<1000.
local p,L;
L:=[];
for p in Primes do
 if Int((p+2)/7)=((p+2)/7) and p^3<n then
   L:=Concatenation(L,[p^3]);
 fi;
 if Int((p-2)/7)=((p-2)/7) and p^3<n then
   L:=Concatenation(L,[p^3]);
 fi;
 if Int((p+3)/7)=((p+3)/7) and p^3<n then
   L:=Concatenation(L,[p^3]);
 fi;
 if Int((p-3)/7)=((p-3)/7) and p^3<n then
   L:=Concatenation(L,[p^3]);
 fi;
od;
return L;
end;

ram_module_hurwitz:=function(p)
##
## input: p is a Hurwitz prime; output: [m1,...,mn],
## where n = # conj classes of G=PSL(2,p)
## and mi = mult of pi_i in ram mod of Hurwitz crv with
## aut gp G.
## Here Irr(G) = [pi_1,...,pi_n] (in that order).
##
local A,B,C,D,pi,G,i,j,n,n0,H,G1,H_chars,CG,G_chars,w,m,theta,pi_theta_to_the;
  if not(p in hurwitz_primes(500)) then
      Print("Input is not a (small) Hurwitz prime.\n");
      return [];
  fi;
  G:=PSL(2,p);
  H:=[];
  H_chars:=[];
  n0:=[ [ Z(p)^0, Z(p) ], [ Z(p)*0, Z(p)^0 ] ];
  CG:=ConjugacyClassesSubgroups(G);
  H[1]:=Representative(CG[2]); # size 2
  H[2]:=Representative(CG[3]); # size 3
  n :=Size(CG);
  for i in [1..n] do
    if Size(Representative(CG[i]))=7 then
      H[3]:=Representative(CG[i]);
    fi;
  od;
## H[3]:=Group(n0);       # size 7
  for i in [1..Size(H)] do
   H_chars[i]:=Irr(H[i]);
  od;
  G_chars:=Irr(G);
  m:=[];
  m[1]:=[];m[2]:=[];m[3]:=[];
  theta:=List([1..3],i->H_chars[i][2]);
  pi_theta_to_the:=List([1..3],i->Sum([1..(Size(H[i])-1)],
      j->j*InducedClassFunction(theta[i]^(j),G)));
#Print("\n\n pi_theta_to_the = ",Sum(pi_theta_to_the),"\n\n");
  for i in [1..3] do
    m[i]:=List(G_chars, pi->ScalarProduct(pi_theta_to_the[i],pi))/Size(H[i]);
  od;
#Print("\n\n m = ",m,"\n\n");
#Print("Multiplicities (by definition):\n ",Sum(m),"\n");
return Sum(m);
####
n:=Length(G_chars);
A:=[];
B:=[];
C:=[];
D:=[];
for i in [1..n] do
  pi:=G_chars[i];
  A[i]:=ScalarProduct(H_chars[1][1],RestrictedClassFunction(pi,H[1]));
  B[i]:=ScalarProduct(H_chars[2][1],RestrictedClassFunction(pi,H[2]));
  C[i]:=ScalarProduct(H_chars[3][1],RestrictedClassFunction(pi,H[3]));
  D[i]:=DegreeOfCharacter(pi);
od;
#Print("Multiplicities (predicted by JK - assumes all R_\ell=1):\n ",(1/2)*(3*D-A-B-C),"\n");
end;

#########example
#gap> ram_module_hurwitz(13); time;
#
#m = [ [0,2,2,3,3,3,3,4,3],[0,2,2,4,4,4,4,5,5],[0,3,3,5,5,5,6,6,6] ]
#
#Multiplicities (by definition):
# [ 0, 7, 7, 12, 12, 12, 13, 15, 14 ]
#Multiplicities (predicted by JK - assumes all R_ell=1):
# [ 0, 7, 7, 12, 12, 12, 13, 15, 14 ]
#1640
#

ram_module_X:=function(p)
## p is a Hurwitz prime
## output is (m1,...,mn)
## where n = # conj classes of G=PSL(2,p)
## and mi = mult of pi_i in ram mod of G
##
local G,i,j,n,n0,H,G1,H_chars,CG,G_chars,w,m,theta,pi_theta_to_the;
  if not(p in hurwitz_primes(500)) then
      Print("Input is not a (small) Hurwitz prime.\n");
      return [];
  fi;
  G:=PSL(2,p);
  H:=[];
  H_chars:=[];
  n0:=[ [ Z(p)^0, Z(p) ], [ Z(p)*0, Z(p)^0 ] ];
  CG:=ConjugacyClassesSubgroups(G);
  H[1]:=Representative(CG[2]); # size 2
  H[2]:=Representative(CG[3]); # size 3
  n :=Size(CG);
  for i in [1..n] do
    if Size(Representative(CG[i]))=7 then
      H[3]:=Representative(CG[i]);
    fi;
  od;
## H[3]:=Group(n0);       # size 7
  for i in [1..Size(H)] do
   H_chars[i]:=Irr(H[i]);
  od;
  G_chars:=Irr(G);
  m:=[];
  m[1]:=[];m[2]:=[];m[3]:=[];
  theta:=List([1..3],i->H_chars[i][2]);
  pi_theta_to_the:=List([1..3],i->Sum([1..(Size(H[i])-1)],
      j->j*InducedClassFunction(theta[i]^(j),G)));
#Print("\n\n pi_theta_to_the = ",Sum(pi_theta_to_the),"\n\n");
  for i in [1..3] do
    m[i]:=List(G_chars, pi->ScalarProduct(pi_theta_to_the[i],pi))/Size(H[i]);
  od;
#Print("\n\n m = ",m,"\n\n");
return Sum(m);
end;

ram_module_X_JK:=function(p)
## p is a Hurwitz prime
## output is (m1,...,mn)
## where n = # conj classes of G=PSL(2,p)
## and mi = "mult of pi_i in ram mod of G" using JK formula
##
local G,i,j,n,n0,H,G1,H_chars,CG,G_chars,A,B,C,D,pi;
G:=PSL(2,p);
H:=[];
H_chars:=[];
n0:=[ [ Z(p)^0, Z(p) ], [ Z(p)*0, Z(p)^0 ] ];
CG:=ConjugacyClassesSubgroups(G);
H[1]:=Representative(CG[2]); # size 2
H[2]:=Representative(CG[3]); # size 3
n:=Size(CG);
for i in [1..n] do
  if Size(Representative(CG[i]))=7 then
    H[3]:=Representative(CG[i]);
  fi;
od;
## H[3]:=Group(n0);       # size 7
for i in [1..Size(H)] do
 H_chars[i]:=Irr(H[i]);
od;
G_chars:=Irr(G);
n:=Length(G_chars);
A:=[];
B:=[];
C:=[];
D:=[];
for i in [1..n] do
  pi:=G_chars[i];
  A[i]:=ScalarProduct(H_chars[1][1],RestrictedClassFunction(pi,H[1]));
  B[i]:=ScalarProduct(H_chars[2][1],RestrictedClassFunction(pi,H[2]));
  C[i]:=ScalarProduct(H_chars[3][1],RestrictedClassFunction(pi,H[3]));
  D[i]:=DegreeOfCharacter(pi);
od;
return (1/2)*(3*D-A-B-C); #assumes all R_\ell=1
end;

equiv_deg_module_hurwitz:=function(p,ii,r)
## p is a Hurwitz prime
## ii  = 1 for H[1], size 2
## ii  = 2 for H[2], size 3
## ii  = 3 for H[3], size p
## output is (m1,...,mn)
## where n = # conj classes of G=PSL(2,p)
## and mi = mult of pi_i in deg_equiv module of G
##
local G,i,j,n,n0,H,G1,H_chars,CG,G_chars,w,m,theta,pi_theta,pi_theta_to_the;
G:=PSL(2,p);
H:=[]; # 3 decomp gps
H_chars:=[];
n0:=[ [ Z(p)^0, Z(p) ], [ Z(p)*0, Z(p)^0 ] ];
CG:=ConjugacyClassesSubgroups(G);
H[1]:=Representative(CG[2]); # size 2
H[2]:=Representative(CG[3]); # size 3
n:=Size(CG);
for i in [1..n] do
  if Size(Representative(CG[i]))=7 then
    H[3]:=Representative(CG[i]);
  fi;
od;
## H[3]:=Group(n0);       # size 7
for i in [1..3] do
 H_chars[i]:=Irr(H[i]);
od;
G_chars:=Irr(G);
m:=[];
m[1]:=[];m[2]:=[];m[3]:=[];
theta:=List([1..3],i->H_chars[i][2]);
pi_theta_to_the:=List([1..3],i->List([1..r],j->InducedClassFunction(theta[i]^(-j),G)));
for i in [1..3] do
 for j in [1..r] do
  m[i][j]:=List(G_chars, pi->ScalarProduct(pi_theta_to_the[i][j],pi));
 od;
od;
return Sum(m[ii]);
end;

# equiv_deg_module_X(7,2,1);
#[ 0, 1, 1, 2, 2, 3 ]
# equiv_deg_module_X(7,1,1);
#[ 0, 2, 2, 2, 4, 4 ]
# equiv_deg_module_X(7,3,1);
#[ 0, 1, 0, 1, 1, 1 ]
# equiv_deg_module_X(7,3,7);
#[ 1, 3, 3, 6, 7, 8 ]
# equiv_deg_module_X(7,1,2);
#[ 1, 3, 3, 6, 7, 8 ]

cuberoots:=function(q)
local x,L;
 L:=[];
 for x in GF(q^2) do
   if x<>Zero(GF(q^2)) and x<>One(GF(q^2)) and x^3=One(GF(q^2)) then
     L:=Concatenation(L,[x]);
   fi;
 od;
 return L;
end;

#gap> cuberoots(41);
#[ Z(41^2)^560, Z(41^2)^1120 ]
#gap> cuberoots(43);
#[ Z(43)^14, Z(43)^28 ]
#gap> fourthroots(43);
#[ Z(43^2)^462, Z(43^2)^1386 ]
#gap> fourthroots(41);
#[ Z(41)^10, Z(41)^30 ]
#gap> fourteenthroots(41);
#[ Z(41^2)^120, Z(41^2)^600, Z(41^2)^1560, Z(41^2)^1320, Z(41^2)^1080, Z(41^2)^360 ]
#gap> fourteenthroots(43);
#[ Z(43)^3, Z(43)^9, Z(43)^15, Z(43)^27, Z(43)^33, Z(43)^39 ]

fourteenthroots:=function(q)
local x,L;
 L:=[];
 for x in GF(q^2) do
   if x<>Zero(GF(q^2)) and x<>One(GF(q^2)) and x^2<>One(GF(q^2))
                     and x^7<>One(GF(q^2)) and x^(14)=One(GF(q^2)) then
     L:=Concatenation(L,[x]);
   fi;
 od;
 return L;
end;

fourthroots:=function(q)
local x,L;
 L:=[];
 for x in GF(q^2) do
   if x<>Zero(GF(q^2)) and x<>One(GF(q^2)) and
      x^2<>One(GF(q^2)) and x^(4)=One(GF(q^2)) then
     L:=Concatenation(L,[x]);
   fi;
 od;
 return L;
end;

Nalpha:=function(alpha,q)
 # alpha a character of GF(q)^x
 local root3,root4,root14,count;
 count := 0;
 root3:=cuberoots(q)[1];
 root4:=fourthroots(q)[1];
 root14:=fourteenthroots(q)[1];
 if root3 in GF(q) and root3^alpha<>1 then count:=count+1; fi;
 if root4 in GF(q) and root4^alpha<>1 then count:=count+1; fi;
 if root14 in GF(q) and root14^alpha<>1 then count:=count+1; fi;
 return count;
end;

Nbeta:=function(beta,q)
 # beta a character of GF(q^2)^x
 local root3,root4,root14,count;
 count := 0;
 root3:=cuberoots(q)[1];
 root4:=fourthroots(q)[1];
 root14:=fourteenthroots(q)[1];
 if root3 in GF(q) and root3^beta<>1 then count:=count+1; fi;
 if root4 in GF(q) and root4^beta<>1 then count:=count+1; fi;
 if root14 in GF(q) and root14^beta<>1 then count:=count+1; fi;
 return count;
end;

ext: gap

max_stars_repo_name: bopopescu/sage

max_issues_repo_name: Ivo-Maffei/sage

max_forks_repo_name: dimpase/sage


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

Class(CScratchUnparserProg, CUnparser, rec(
    dma_signal := (self,o,i,is) >> Print(Blanks(i), self.opts.dmaSignal(self.opts), self.pinfix(o.args, ", "), ";\n"),
    dma_wait := (self,o,i,is) >> Print(Blanks(i), self.opts.dmaWait(self.opts), "", self.pinfix(o.args, ", "), ";\n"),

    cpu_signal := (self,o,i,is) >> Print(Blanks(i), self.opts.cpuSignal(self.opts), self.pinfix(o.args, ", "), ";\n"),
    cpu_wait := (self,o,i,is) >> Print(Blanks(i), self.opts.cpuWait(self.opts), "", self.pinfix(o.args, ", "), ";\n"),

    dma_fence := (self,o,i,is) >> Print(Blanks(i), self.opts.dmaFence(self.opts), "();\n"),

    dma_load := (self,o,i,is) >> Print(Blanks(i), self.opts.dmaLoad(self.opts),
        "(", self(o.loc,i,is), ", ", self(o.exp,i,is), ", ", self(o.size,i,is), ");\n"),

    dma_store := (self,o,i,is) >> Print(Blanks(i), self.opts.dmaStore(self.opts),
        "(", self(o.loc,i,is), ", ", self(o.exp,i,is), ", ", self(o.size,i,is), ");\n"),

    par_exec := (self,o,i,is) >> Print(Blanks(i), "parallel {\n", DoForAll(o.cmds, c -> self(c, i+is, is)), Blanks(i), "}\n"),

    decl := meth(self,o,i,is)
        local arrays, memarrays, scratcharrays, romarrays, other, l, arri, myMem;
        [arrays, other] := SplitBy(o.vars, x->IsArray(x.t));
        [memarrays, arrays] := SplitBy(arrays, i->IsBound(i.t.qualifiers) and self.opts.memModifier in i.t.qualifiers);
        [scratcharrays, arrays] := SplitBy(arrays, i->IsBound(i.t.qualifiers) and self.opts.scratchModifier in i.t.qualifiers);
        [romarrays, arrays] := SplitBy(arrays, i->IsBound(i.t.qualifiers) and self.opts.romModifier in i.t.qualifiers);

        if Length(arrays) > 0 then
            DoForAll(arrays, v -> Print(Blanks(i), self.opts.arrayBufModifier, " ", self.declare(v.t, v, i, is), ";\n"));
        fi;
        if Length(memarrays) > 0 then
            DoForAll(memarrays, v -> Print(Blanks(i), self.opts.arrayBufModifier, " ", self.opts.memModifier, " ", self.declare(v.t, v, i, is), ";\n"));
        fi;
        if Length(scratcharrays) > 0 then
            DoForAll(scratcharrays, v -> Print(Blanks(i), self.opts.arrayBufModifier, " ", self.opts.scratchModifier, " ", self.declare(v.t, v, i, is),";\n"));
        fi;
        if Length(romarrays) > 0 then
            DoForAll(romarrays, v -> Print(Blanks(i), self.opts.arrayBufModifier, " ", self.opts.romModifier, " ", self.declare(v.t, v, i, is), ";\n"));
        fi;

        if (Length(other)>0) then
            other:=SortRecordList(other,x->x.t);
            for l in other do
               Sort(l, (a,b)->a.id < b.id);
               Print(Blanks(i), self.declare(l[1].t, l, i, is), ";\n");
            od;
        fi;

        self(o.cmd, i, is);

        #Pop arena for this decl
        if IsBound(self.opts.useMemoryArena) and self.opts.useMemoryArena and Length(arrays) > 0 and arrays[1].id[1] <> 'D' then
          myMem := 0;
          for arri in arrays do
             # Account for vector allocations in memory arena (which is scalar)
             myMem := myMem + (arri.t.size * When(IsBound(arri.t.t) and ObjId(arri.t.t)=TVect, arri.t.t.size, 1));
          od;
          if ObjId(myMem) = Value then myMem := myMem.v; fi;
          Print(Blanks(i));
          Print("arenalevel += ", myMem, ";\n" );
        fi;
    end,

    swp_loop := (self, o, i, is) >> let(v := o.var, lo := o.range[1], hi := Last(o.range),
        Print(When(IsBound(self.opts.looppragma), self.opts.looppragma(o,i,is)),
          Blanks(i), "for(int ", v, " = ", lo, "; ", v, " <= ", hi, "; ", v, "++) { // SWP loop\n",
          self(o.cmd,i+is,is),
          Blanks(i), "}\n")),
));


text:
# Copyright (c) 2018-2020, Carnegie Mellon University
# See LICENSE for details
#
# Contains the rewrite rulesets for SPL-expression simplification
# After the Breakdown Stage
# Current Demo version focuses on CNOTs, and so we do not cancel single-qubit gates
# Eventually we will convert all gates to generic rotations in Quantum_Format
# Combine rotations in Quantum_Simplify, and convert back into non-generic rotations in Quantum_Terminate
#
# Additionally, Junction steps should be factored out in this stage versus in the unparser to enable a
# few more optimizations. A SWE effort is underway.
#

##
#F CheckIdentitySplit( <a>, <b> )
##
## return true if we can split an Identity transform into smaller matrices
CheckIdentitySplit := function(a, b)
    if(ObjId(a) = I) then
        if a.params[1] > 2 then
            return true;
        fi;
    fi;
    if(ObjId(b) = I) then
        if b.params[1] > 2 then
            return true;
        fi;
    fi;
    return false;
end;

##
#F CheckIdentity( <a>, <b> )
##
## are both objects identity
CheckIdentity := function(a, b)
    if(ObjId(a) = I) then
        if(ObjId(b) = I) then
            return true;
        fi;
    fi;
    return false;
end;

##
#F PullInTensor( <ch1>, <ch2> )
##
## compress two tensor expressions by pulling in symbols wherever sizes align
PullInTensor := function(ch1, ch2)
    local pos1, pos2, e1, e2, up1, up2, siz1, siz2;
    pos1 := 1;
    pos2 := 1;
    siz1 := 0;
    siz2 := 0;
    while pos1 <= Length(ch1) and pos2 <= Length(ch2) do
        e1 := ch1[pos1];
        e2 := ch2[pos2];
        if ( (e1.dims()[1] = e2.dims()[1]) and (pos1 = pos2) ) then
            ch1[pos1] := ch1[pos1] * ch2[pos2];
            ch2[pos2] := I(e2.dims()[1]);
        fi;
        up1 := siz1 + log(e1.dims()[1], 2).v;
        up2 := siz2 + log(e2.dims()[1], 2).v;
        if (up1 <= up2) then
            pos1 := pos1 + 1;
            siz1 := siz1 + log(e1.dims()[1], 2).v;
        fi;
        if (up1 > up2) then
            pos2 := pos2 + 1;
            siz2 := siz2 + log(e2.dims()[1], 2).v;
        fi;
    od;
    return [ch1, ch2];
end;

##
#F CanPullIn( <ch1>, <ch2> )
##
## Similar to PullInTensor, but just checks to see if there is any simplification possible
CanPullIn := function(ch1, ch2)
    local pos1, pos2, e1, e2, up1, up2, siz1, siz2;
    pos1 := 1;
    pos2 := 1;
    siz1 := 0;
    siz2 := 0;
    while pos1 <= Length(ch1) and pos2 <= Length(ch2) do
        e1 := ch1[pos1];
        e2 := ch2[pos2];
        if ( (e1.dims()[1] = e2.dims()[1]) and (pos1 = pos2) and ObjId(e2) <> I) then
            return true;
        fi;
        up1 := siz1 + log(e1.dims()[1], 2).v;
        up2 := siz2 + log(e2.dims()[1], 2).v;
        if (up1 <= up2) then
            pos1 := pos1 + 1;
            siz1 := siz1 + log(e1.dims()[1], 2).v;
        fi;
        if (up1 > up2) then
            pos2 := pos2 + 1;
            siz2 := siz2 + log(e2.dims()[1], 2).v;
        fi;
    od;
    return false;
end;

##
#F ValidTensSimplifyMatch( <e1>, <e2> )
##
## SWrapper function to deterrmine if the tensor rewrite rule can be applied
ValidTensSimplifyMatch := function(e1, e2)
    if ObjId(e1)=Tensor and ObjId(e2)=Tensor then
        if CanPullIn(e1._children, e2._children) then
            return true;
        fi;
    fi;
    return false;
end;

##
#F CombineReord( <ord1>, <dir1>, <ord2>, <dir2>, <arch> )
##
## Combine Reorder steps
CombineReord := function(ord1, dir1, ord2, dir2, arch)
    local idx, adj1, adj2, adj3, start, n, t1, t2, a, b, new_list, o1, o2;
    adj1 := [];
    adj2 := [];
    if dir1 = -1 then
        idx := 0;
        for o1 in ord1 do
            Add(adj1, [idx, o1]);
            idx := idx + 1;
        od;
    fi;
    if dir1 = 1 then
        idx := 0;
        for o1 in ord1 do
            Add(adj1, [o1, idx]);
            idx := idx + 1;
        od;
    fi;
    if dir2 = -1 then
        idx := 0;
        for o2 in ord2 do
            Add(adj2, [idx, o2]);
            idx := idx + 1;
        od;
    fi;
    if dir2 = 1 then
        idx := 0;
        for o2 in ord2 do
            Add(adj2, [o2, idx]);
            idx := idx + 1;
        od;
    fi;
    # now have adj1 and adj2
    adj3 := [];
    for a in adj1 do
        start := a[1];
        t1 := a[2];
        t2 := -1;
        for b in adj2 do
            if b[1] = t1 then
                t2 := b[2];
            fi;
        od;
        Add(adj3, [start, t2]);
    od;
    # full combine, and then check validity. If it isnt valid, then just dont simplify
    # TODO do a partial recombine if the entire reorder cannot be combined and be still implementable
    # not much harder, just a lot of code
    new_list := [0..Length(ord1)-1];
    idx := 1;
    for n in new_list do
        for a in adj3 do
            if n = a[1] then
                new_list[idx] := a[2];
            fi;
        od;
        idx := idx + 1;
    od;
    if VerifyPath(new_list, [0..Length(ord1)-1], arch) = true then
        return Reord(new_list, arch, 1);
    fi;
    return (Reord(ord1, arch, dir1) * Reord(ord2, arch, dir2));
end;

##
#F Quantum_Format ruleset
##
Class(Quantum_Format, RuleSet);
RewriteRules(Quantum_Format, rec(
    # Remove Grp structure
    remove_grp := ARule( Grp, [@(1)], x-> [(@(1).val)] ),

    # Flatten Tensors
    flatten_tensor_asoc := ARule(Tensor, [@(1), @(2).cond(e-> ObjId(@(1).val)=Tensor or ObjId(e)=Tensor)],
	    e -> let(
	           ch1 := Cond(ObjId(@(1).val)=Tensor, @(1).val._children, [@(1).val]),
	           ch2 := Cond(ObjId(@(2).val)=Tensor, @(2).val._children, [@(2).val]),
	           ch1::ch2
	         )),

    #Split identities into I(2)
    split_identity := ARule(Tensor, [@(1),  @(2).cond(e-> (CheckIdentitySplit(@(1).val, e)) )],
	    e -> let(
	           ch1 := Cond(ObjId(@(1).val)=I and @(1).val.params[1] > 2, Tensor(I(2), I(@(1).val.params[1]/2)), [@(1).val]),
	           ch2 := Cond(ObjId(@(2).val)=I and @(2).val.params[1] > 2, Tensor(I(2), I(@(2).val.params[1]/2)), [@(2).val]),
               [ch1]::[ch2]
	         )),
));

##
#F Quantum_Simplify ruleset
##
Class(Quantum_Simplify, RuleSet);
RewriteRules(Quantum_Simplify, rec(

    # Pull-In Tensors
    pull_in_tensor := ARule(Compose, [@(1), @(2).cond(e-> ValidTensSimplifyMatch(@(1).val, e))],
	    e -> let(
               # both 1 ans 2 are tensors, can access a list of arguments via _children
               ch := PullInTensor(@(1).val._children, @(2).val._children),
	           [Tensor(ch[1])]::[Tensor(ch[2])]
	         )),

    # Remove Identity Multiplications
    remove_identity_right := ARule( Compose, [@(1), @(2).cond(e -> ObjId(e)=I)], x-> [(@(1).val)]),
    remove_identity_left := ARule( Compose, [@(1).cond(e -> ObjId(e)=I), @(2)], x-> [(@(2).val)]),

    # Cancel Composition of CNOTS
    cancel_cnot := ARule(Compose, [[CNOT, @(1), @(2)], [CNOT, @(3).cond(x -> (x = @(1).val)), @(4).cond(x -> (x = @(2).val))]], x->[I(2^(@(1).val + 1))] ),

    # Combine Reorder Compositions
    cancel_reorder := ARule(Compose, [[Reord, @(1), @(2), @(3)], [Reord, @(4), @(5), @(6)]], x-> [CombineReord(@(1).val, @(3).val, @(4).val, @(6).val, @(5).val)]),
));

##
#F Quantum_Terminate ruleset
##
Class(Quantum_Terminate, RuleSet);
RewriteRules(Quantum_Terminate, rec(
    # Recombine Identities - Good
    recomb_identity := ARule(Tensor, [@(1),  @(2).cond(e-> (CheckIdentity(@(1).val, e) ))],
	    e -> [I(@(1).val.params[1] * @(2).val.params[1])]),

    # Remove Identity Multiplications - Good
    remove_identity_right := ARule( Compose, [@(1), @(2).cond(e -> ObjId(e)=I)], x-> [(@(1).val)]),
    remove_identity_left := ARule( Compose, [@(1).cond(e -> ObjId(e)=I), @(2)], x-> [(@(2).val)]),

));

##
#F QuantumRewriteInternal( <s>, <opts> )
##
## Rewrites a quantum SPL expression, for internal use
QuantumRewriteInternal := function (s, opts)
    s := ApplyStrategy(s, [Quantum_Format], BUA, opts);
    s := ApplyStrategy(s, [Quantum_Simplify], BUA, opts);
    s := ApplyStrategy(s, [Quantum_Terminate], BUA, opts);
    return s;
end;

##
#F BestCircuitInternal( <t>, <opts> )
##
## Finds the best circuit implementing transform t, for internal use
BestCircuitInternal := function (t, opts)
    local dpopts, best, s;
    dpopts := rec(verbosity := 0, hashTable := HashTableDP());
    dpopts.measureFunction := (rt, opts) ->
        let(c := SPLRuleTree(rt), # generate SPL Ruletree and then run a collect on CNOTs
            c2 := QuantumRewriteInternal(c, opts),
            Length(Collect(c2, @(1,CNOT)))
            );
    best := DP(t, dpopts, opts);
    s := QuantumRewriteInternal(SPLRuleTree(best[1].ruletree), opts);
    return s;
end;

##
#F TerminateReord( <l>, <arch>, <dir> )
##
## Convert a reorder object to CNOTs
TerminateReord := function (l, arch, dir)
    local mats, mat, dir, spl;
    mats := ShiftMat(l, [0..Length(arch)-1], arch, Length(arch));
    if dir = -1 then
        mat := mats[2];
    fi;
    if dir = 1 then
        mat := mats[1];
    fi;
    # now mat in an unreduced Non-terminal expression
    # we need the SPL expression
    spl := BestCircuitInternal(mat, SpiralDefaults);
    return spl;
end;

##
#F Quantum_Reorder ruleset
##
Class(Quantum_Reorder, RuleSet);
RewriteRules(Quantum_Reorder, rec(

    # Combine Reorder Compositions
    reorder_convert := Rule([Reord, @(1), @(2), @(3)], x-> TerminateReord(@(1).val, @(2).val, @(3).val)),
));

# ------------------------- For External Use ----------------------------------------

##
#F QuantumRewrite( <s>, <opts> )
##
## Simplifies an SPL expression representing a quantum circuit (s)
QuantumRewrite := function (s, opts)
    s := ApplyStrategy(s, [Quantum_Format], BUA, opts);
    s := ApplyStrategy(s, [Quantum_Simplify], BUA, opts);
    s := ApplyStrategy(s, [Quantum_Reorder], BUA, opts);
    s := RulesSums(s);
    s := ApplyStrategy(s, [Quantum_Format], BUA, opts);
    s := ApplyStrategy(s, [Quantum_Simplify], BUA, opts);
    s := ApplyStrategy(s, [Quantum_Terminate], BUA, opts);
    return s;
end;

##
#F BestCircuit( <t>, <opts> )
##
## Finds the best circuit implementing transform t
## Hardcoded to count number CNOT gates as the cost function
BestCircuit := function (t, opts)
    local dpopts, best, s, cnot_list;
    dpopts := rec(verbosity := 1, hashTable := HashTableDP());
    dpopts.measureFunction := (rt, opts) ->
        let(c := SPLRuleTree(rt), # generate SPL Ruletree and then run a collect on CNOTs
            c2 := QuantumRewrite(c, opts),
            Length(Collect(c2, @(1,CNOT)))
            );
    best := DP(t, dpopts, opts);
    Print("\n SPL \n");
    Print(SPLRuleTree(best[1].ruletree));
    Print("\n");
    s := QuantumRewrite(SPLRuleTree(best[1].ruletree), opts);
    Print("\ncost :");
    Print(Length(Collect(s, @(1,CNOT))));
    Print("\n");
    Print("---------------- BestCircuit: \n");
    return s;
end;

max_stars_repo_name: spiral-software/spiral-package-quantum

max_issues_repo_name: spiral-software/spiral-package-quantum

max_forks_repo_name: spiral-software/spiral-package-quantum


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

Import(platforms, profiler, code);
Import(platforms.sse, platforms.avx, platforms.intel);

Class(OSInfo, rec(osname := "",
    is8Bit := False,
    is16bit := False,
    is32bit := False,
    is64bit := False,
    isWindows := False,
    isLinux := False,
    isCygwin := False,
    isDarwin := False,
    info := self >> Cond(IsBound(self.vendor), Print(self.vendor, " ", self.osname), Print(self.osname)),
    useColor := True
));

Class(OSWindows, OSInfo, rec(
    setTitle := meth(self, title)
        Exec(Concat("title ", title));
    end,
    vendor := "Microsoft",
    isWindows := True,
    useColor := False
));

Class(OSWindows32, OSWindows, rec(is16bit := True, is32bit := True));
Class(OSWindows64, OSWindows, rec(is16bit := True, is32bit := True, is64bit :=True));
Class(OSLinux32, OSInfo, rec(osname := "Linux32", is32bit := True, isLinux := True));
Class(OSLinux64, OSInfo, rec(osname := "Linux64", is32bit := True, is64bit := True, isLinux := True));
Class(OSArmLinux, OSInfo, rec(osname := "GNU/Linux", is16bit := True, is32bit := True, isLinux := True));
Class(OSCygwin32, OSInfo, rec(osname := "Cygwin32", is32bit := True, isCygwin := True));
Class(OSDarwin, OSInfo, rec(osname := "OSX/Darwin", is32bit := True, is64bit := True, isDarwin := True));

SupportedOSs := rec(
    WindowsNT4 := CopyFields(OSWindows32, rec(osname := "Windows NT 4.0")),
    Windows2000 := CopyFields(OSWindows32, rec(osname := "Windows 2000")),
    WindowsXP32 := CopyFields(OSWindows32, rec(osname := "WindowsXP 32-bit")),
    WindowsXP64 := CopyFields(OSWindows64, rec(osname := "WindowsXP 64-bit")),
    WindowsVista := CopyFields(OSWindows64, rec(osname := "Windows Vista")),
    Windows7 := CopyFields(OSWindows64, rec(osname := "Windows 7")),
    Windows8 := CopyFields(OSWindows64, rec(osname := "Windows 8")),
    Linux32 := OSLinux32,
    Linux64 := OSLinux64,
    Cygwin32 := OSCygwin32,
    Darwin := OSDarwin,
    ArmLinux := OSArmLinux,
);

Class(CPUInfo, rec(
    hasDouble := True,
    hasFloat := True,
    intSize := 32,
    cores := 1,
    cpuname := "",
    vendor := "",
    info := self >> Chain(
        Print(Cond(IsBound(self.vendor), Print(self.vendor, " ", self.cpuname), Print(self.osname))),
        When(IsBound(self.freq) and self.freq > 0, Print(" at ", self.freq, " MHz")),
        When(self.cores > 1, Print(", ", self.cores, " cores")),
        When(IsBound(self.SIMDname), Print(", ", self.SIMDname))
    )
));

Class(IntelCPU, CPUInfo, rec(
    vendor := "Intel",
    getSimdIsa := (self, dt) >> When(self.SIMD().hasAVX(),
        Cond(
            dt = T_Real(32), AVX_8x32f,
            dt = T_Real(64), AVX_4x64f,
            dt),
        Cond(
            dt = T_Real(32), SSE_4x32f,
            dt = T_Real(64), SSE_2x64f,
            dt)
    ),
    getOpts := arg >> IAGlobals.getOpts(Drop(arg, 1))
));

Class(AMDCPU, CPUInfo, rec(
    vendor := "AMD"
));

Class(STICPU, CPUInfo, rec(
    vendor := "STI"
));

Class(PowerPC, CPUInfo, rec(
    vendor := "IBM"
));

Class(DPA, CPUInfo, rec(
       vendor := "CMU"
));

Class(ARM, CPUInfo, rec(
       vendor := "Raspberry Pi Foundation"
));

SupportedCPUs := rec(
    Pentium := rec(),
    PentiumPro := rec(),
    PentiumII := rec(),
    PentiumIII := rec(),
    Pentium4 := Class(Pentium4, IntelCPU, rec(
        cpuname := "Pentium4",
        SIMD := () -> CopyFields(platforms.SIMDArchitectures, rec(hasMMX := True, hasSSE:=True, hasSSE2:=True, hasSSE3:=False)),
        SIMDname := "SSE2",
        cores := 1,
        default_lang := "c.icl.opt.pentium4"
    )),
    Pentium4Extreme := Class(Pentium4Extreme, IntelCPU, rec(
        cpuname := "Pentium4Extreme",
        SIMD := () -> CopyFields(platforms.SIMDArchitectures, rec(hasMMX := True, hasSSE:=True, hasSSE2:=True, hasSSE3:=True)),
        SIMDname := "SSE3",
        cores := 1,
        default_lang := "c.icl.opt.pentium4extreme"
    )),
    PentiumD := rec(),
    Xeon := Class(Xeon, IntelCPU, rec(
        cpuname := "Xeon",
        SIMD := () -> CopyFields(platforms.SIMDArchitectures, rec(hasMMX := True, hasSSE:=True, hasSSE2:=True, hasSSE3:=True, hasSSSE3:=True, hasSSE4_1 := True, hasSSE4_2 := True)),
        SIMDname := "SSE4.2",
        cores := 4,
        profile := rec(
            EM64T := (arg) -> CopyFields(default_profiles.win_x64_icc, rec(CFLAGS := arg -> "/O3 /G7 /QxS")),
            IA32 := (arg) -> CopyFields(default_profiles.win_x86_icc, rec(CFLAGS := arg -> "/O3 /G7 /QxS")),
            threads := (arg) -> CopyFields(default_profiles.win_x86_icc_threads, rec(CFLAGS := arg -> "/O3 /G7 /QxS"))
        ),
        default_lang := "c.icl.opt.corei7",
        smp_lang := "c.icl.smp_corei7",
        OpenMP_lang := "c.icl.openmp_corei7"
    )),
    XeonMP := rec(),
    PentiumM := Class(PentiumM, IntelCPU, rec(
        cpuname := "PentiumM",
        SIMD := () -> CopyFields(platforms.SIMDArchitectures, rec(hasMMX := True, hasSSE:=True, hasSSE2:=True, hasSSE3:=False)),
        SIMDname := "SSE2",
        cores := 1,
        default_lang := "c.icl.opt.pentiumM"
    )),
    CoreDuo := Class(CoreDuo, IntelCPU, rec(
        cpuname := "CoreDuo",
        SIMD := () -> CopyFields(platforms.SIMDArchitectures, rec(hasMMX := True, hasSSE:=True, hasSSE2:=True, hasSSE3:=True)),
        SIMDname := "SSE3",
        cores := 2,
        profile := rec(
            IA32 := (arg) -> CopyFields(default_profiles.win_x86_icc, rec(CFLAGS := arg -> "/O3 /G7 /QxKWP")),
            threads := (arg) -> When(LocalConfig.osinfo.isWindows(),
        CopyFields(default_profiles.win_x86_icc_threads, rec(CFLAGS := arg -> "/O3 /G7 /QxKWP")),
        CopyFields(default_profiles.linux_x86_threads, rec(CFLAGS := arg -> "-msse3"))
)
        ),
        default_lang := "c.icl.opt.core",
        smp_lang := "c.icl.smp_core",
        OpenMP_lang := "c.icl.openmp_core"
    )),
    Core2Duo := Class(Core2Duo, IntelCPU, rec(
        cpuname := "Core2Duo",
        SIMD := () -> CopyFields(platforms.SIMDArchitectures, rec(hasMMX := True, hasSSE:=True, hasSSE2:=True, hasSSE3:=True, hasSSSE3:=True)),
        SIMDname := "SSSE3",
        cores := 2,
        profile := rec(
            EM64T := (arg) -> CopyFields(default_profiles.win_x64_icc, rec(CFLAGS := arg -> "/O3 /G7 /QxSSSE3")),
            IA32 := (arg) -> CopyFields(default_profiles.win_x86_icc, rec(CFLAGS := arg -> "/O3 /G7 /QxKWP")),
            threads := (arg) -> CopyFields(default_profiles.win_x86_icc_threads, rec(CFLAGS := arg -> "/O3 /G7 /QxKWP"))
        ),
        default_lang := "c.icl.opt.core2",
        smp_lang := "c.icl.smp_core2",
        OpenMP_lang := "c.icl.openmp_core2"
    )),
    Core2Quad := Class(Core2Quad, IntelCPU, rec(
        cpuname := "Core2 Extreme",
        SIMD := () -> CopyFields(platforms.SIMDArchitectures, rec(hasMMX := True, hasSSE:=True, hasSSE2:=True, hasSSE3:=True, hasSSSE3:=True)),
        SIMDname := "SSSE3",
        cores := 4,
        profile := rec(
            EM64T := (arg) -> CopyFields(default_profiles.win_x64_icc, rec(CFLAGS := arg -> "/O3 /G7 /QxSSSE3")),
            IA32 := (arg) -> CopyFields(default_profiles.win_x86_icc, rec(CFLAGS := arg -> "/O3 /G7 /QxKWP")),
            threads := (arg) -> CopyFields(default_profiles.win_x86_icc_threads, rec(CFLAGS := arg -> "/O3 /G7 /QxKWP"))
        ),
        default_lang := "c.icl.opt.core2",
        smp_lang := "c.icl.smp_core2",
        OpenMP_lang := "c.icl.openmp_core2"
    )),

    DPABasic := Class(DPABasic, DPA, rec(
	    cpuname := "DPABasic",
	    cores := 1,
	    SIMD := () -> CopyFields(platforms.SIMDArchitectures, rec( numregs := 32)),
	    )),

    Core2Penryn := Class(Core2Penryn, IntelCPU, rec(
        cpuname := "Core2 Penryn",
        SIMD := () -> CopyFields(platforms.SIMDArchitectures, rec(hasMMX := True, hasSSE:=True, hasSSE2:=True, hasSSE3:=True, hasSSSE3:=True, hasSSE4_1 := True)),
        SIMDname := "SSE4.1",
        cores := 4,
        profile := rec(
            EM64T := (arg) -> CopyFields(default_profiles.win_x64_icc, rec(CFLAGS := arg -> "/O3 /G7 /QxS")),
            IA32 := (arg) -> CopyFields(default_profiles.win_x86_icc, rec(CFLAGS := arg -> "/O3 /G7 /QxS")),
            threads := (arg) -> CopyFields(default_profiles.win_x86_icc_threads, rec(CFLAGS := arg -> "/O3 /G7 /QxS"))
        ),
        default_lang := "c.icl.opt.core2p",
        smp_lang := "c.icl.smp_core2p",
        OpenMP_lang := "c.icl.openmp_core2p"
    )),
    Core_i7 := Class(Core_i7, IntelCPU, rec(
        cpuname := "Core i7",
        SIMD := () -> CopyFields(platforms.SIMDArchitectures, rec(hasMMX := True, hasSSE:=True, hasSSE2:=True, hasSSE3:=True, hasSSSE3:=True, hasSSE4_1 := True, hasSSE4_2 := True)),
        SIMDname := "SSE4.2",
        cores := 4,
        profile := rec(
            EM64T := (arg) -> CopyFields(default_profiles.win_x64_icc, rec(CFLAGS := arg -> "/O3 /G7 /QxS")),
            IA32 := (arg) -> CopyFields(default_profiles.win_x86_icc, rec(CFLAGS := arg -> "/O3 /G7 /QxS")),
            threads := (arg) -> CopyFields(default_profiles.win_x86_icc_threads, rec(CFLAGS := arg -> "/O3 /G7 /QxS"))
        ),
        default_lang := "c.icl.opt.corei7",
        smp_lang := "c.icl.smp_corei7",
        OpenMP_lang := "c.icl.openmp_corei7"
    )),
	# specifically for Spiral FFT GPL 1.0
    Core_AVX := Class(Core_AVX, IntelCPU, rec(
        cpuname := "Intel Core with AVX",
        SIMD := () -> CopyFields(platforms.SIMDArchitectures, rec(hasMMX := True, hasSSE:=True, hasSSE2:=True, hasSSE3:=True, hasSSSE3:=True, hasSSE4_1 := True, hasSSE4_2 := True, hasAVX := True)),
        SIMDname := "AVX",
        cores := 4,
        profile := rec(
            EM64T := (arg) -> CopyFields(default_profiles.win_x64_icc, rec(CFLAGS := arg -> "/O3 /G7 /QxS")),
            IA32 := (arg) -> CopyFields(default_profiles.win_x86_icc, rec(CFLAGS := arg -> "/O3 /G7 /QxS")),
            threads := (arg) -> CopyFields(default_profiles.win_x86_icc_threads, rec(CFLAGS := arg -> "/O3 /G7 /QxS"))
        ),
        default_lang := "c.icl.opt.corei7",
        smp_lang := "c.icl.smp_corei7",
        OpenMP_lang := "c.icl.openmp_corei7"
    )),
	# specifically for Spiral FFT GPL 1.0
    Core_no_AVX := Class(Core_no_AVX, IntelCPU, rec(
        cpuname := "Intel Core without AVX",
        SIMD := () -> CopyFields(platforms.SIMDArchitectures, rec(hasMMX := True, hasSSE:=True, hasSSE2:=True, hasSSE3:=True, hasSSSE3:=True, hasSSE4_1 := True, hasSSE4_2 := True, hasAVX := False)),
        SIMDname := "SSE4.2",
        cores := 4,
        profile := rec(
            EM64T := (arg) -> CopyFields(default_profiles.win_x64_icc, rec(CFLAGS := arg -> "/O3 /G7 /QxS")),
            IA32 := (arg) -> CopyFields(default_profiles.win_x86_icc, rec(CFLAGS := arg -> "/O3 /G7 /QxS")),
            threads := (arg) -> CopyFields(default_profiles.win_x86_icc_threads, rec(CFLAGS := arg -> "/O3 /G7 /QxS"))
        ),
        default_lang := "c.icl.opt.corei7",
        smp_lang := "c.icl.smp_corei7",
        OpenMP_lang := "c.icl.openmp_corei7"
    )),
    Core_i7U := Class(Core_i7U, IntelCPU, rec(
        cpuname := "Core i7 Ultrabook",
        SIMD := () -> CopyFields(platforms.SIMDArchitectures, rec(hasMMX := True, hasSSE:=True, hasSSE2:=True, hasSSE3:=True, hasSSSE3:=True, hasSSE4_1 := True, hasSSE4_2 := True, hasAVX := True)),
        SIMDname := "AVX",
        cores := 2,
        profile := rec(
            EM64T := (arg) -> CopyFields(default_profiles.win_x64_icc, rec(CFLAGS := arg -> "/O3 /G7 /QxS")),
            IA32 := (arg) -> CopyFields(default_profiles.win_x86_icc, rec(CFLAGS := arg -> "/O3 /G7 /QxS")),
            threads := (arg) -> CopyFields(default_profiles.win_x86_icc_threads, rec(CFLAGS := arg -> "/O3 /G7 /QxS"))
        ),
        default_lang := "c.icl.opt.corei7",
        smp_lang := "c.icl.smp_corei7",
        OpenMP_lang := "c.icl.openmp_corei7"
    )),
    Core_i5U := Class(Core_i5U, IntelCPU, rec(
        cpuname := "Core i5U",
        SIMD := () -> CopyFields(platforms.SIMDArchitectures, rec(hasMMX := True, hasSSE:=True, hasSSE2:=True, hasSSE3:=True, hasSSSE3:=True, hasSSE4_1 := True, hasSSE4_2 := True)),
        SIMDname := "SSE4.2",
        cores := 2,
        profile := rec(
            EM64T := (arg) -> CopyFields(default_profiles.win_x64_icc, rec(CFLAGS := arg -> "/O3 /G7 /QxS")),
            IA32 := (arg) -> CopyFields(default_profiles.win_x86_icc, rec(CFLAGS := arg -> "/O3 /G7 /QxS")),
            threads := (arg) -> CopyFields(default_profiles.win_x86_icc_threads, rec(CFLAGS := arg -> "/O3 /G7 /QxS"))
        ),
        default_lang := "c.icl.opt.corei7",
        smp_lang := "c.icl.smp_corei7",
        OpenMP_lang := "c.icl.openmp_corei7"
    )),
    Core_i5_SandyBridge := Class(Core_i5_SandyBridge, IntelCPU, rec(
        cpuname := "Core i5 SandyBridge",
        SIMD := () -> CopyFields(platforms.SIMDArchitectures, rec(hasMMX := True, hasSSE:=True, hasSSE2:=True, hasSSE3:=True, hasSSSE3:=True, hasSSE4_1 := True, hasSSE4_2 := True, hasAVX := True)),
        SIMDname := "AVX",
        cores := 4,
        profile := rec(
            EM64T := (arg) -> CopyFields(default_profiles.win_x64_icc, rec(CFLAGS := arg -> "/O3 /G7 /QxS")),
            IA32 := (arg) -> CopyFields(default_profiles.win_x86_icc, rec(CFLAGS := arg -> "/O3 /G7 /QxS")),
            threads := (arg) -> CopyFields(default_profiles.win_x86_icc_threads, rec(CFLAGS := arg -> "/O3 /G7 /QxS"))
        ),
    )),

    Itanium := rec(),
    Itanium2 := rec(),
    Itanium3 := rec(),
    MPC_G4 := rec(),
    PPC905_G5 := rec(),
    XScale := rec(),
    PowerPC405 := rec(),
    Athlon := rec(),
    AthlonXP := rec(),
    Opteron := rec(),
    OpteronDual := rec(),
    OpteronQuad := Class(OpteronQuad, AMDCPU, rec(
        cpuname := "Dual Opteron 2220",
        SIMD := () -> CopyFields(platforms.SIMDArchitectures, rec(hasMMX := True, hasSSE:=True, hasSSE2:=True, hasSSE3:=True)),
        SIMDname := "SSE3",
        cores := 4,
        profile := rec(
            EM64T := (arg) -> CopyFields(default_profiles.win_x64_icc, rec(CFLAGS := arg -> "/O3 /G7 /QxSSSE3")),
            IA32 := (arg) -> CopyFields(default_profiles.win_x86_icc, rec(CFLAGS := arg -> "/O3 /G7 /QxKWP")),
            threads := (arg) -> CopyFields(default_profiles.win_x86_icc_threads, rec(CFLAGS := arg -> "/O3 /G7 /QxKWP"))
        ),
        default_lang := "c.icl.opt.opteron2220"
    )),
    CellBE := Class(CellBE, STICPU, rec(
        cpuname := "Cell BE",
        SIMD := () -> CopyFields(platforms.SIMDArchitectures, rec(hasSPU := True)),
        SIMDname := "SIMD-SPU",
        cores := 9,
        profile := rec(
            EM64T := (arg) -> CopyFields(default_profiles.win_x64_icc, rec(CFLAGS := arg -> "")),
            IA32 := (arg) -> CopyFields(default_profiles.win_x86_icc, rec(CFLAGS := arg -> "")),
            threads := (arg) -> CopyFields(default_profiles.win_x86_icc_threads, rec(CFLAGS := arg -> ""))
        ),
        default_lang := ""
   )),
   CellBEPS3 := Class(CellBEPS3, STICPU, rec(
        cpuname := "Cell BE (PS3)",
        SIMD := () -> CopyFields(platforms.SIMDArchitectures, rec(hasSPU := True)),
        SIMDname := "SIMD-SPU",
        cores := 6,
        profile := rec(
            EM64T := (arg) -> CopyFields(default_profiles.win_x64_icc, rec(CFLAGS := arg -> "")),
            IA32 := (arg) -> CopyFields(default_profiles.win_x86_icc, rec(CFLAGS := arg -> "")),
            threads := (arg) -> CopyFields(default_profiles.win_x86_icc_threads, rec(CFLAGS := arg -> ""))
        ),
        default_lang := "",
        getSimdIsa := dt -> Cond(
            dt = T_Real(32), platforms.cellSPU.spu_4x32f,
            dt = T_Real(64), platforms.cellSPU.spu_2x64f,
            dt),
   )),
   PowerPC970 := Class(PowerPC970, PowerPC, rec(
        cpuname := "PowerPC 970 (G5)",
        SIMD := () -> CopyFields(platforms.SIMDArchitectures, rec(hasAltiVec := true)),
        SIMDname := "AltiVec",
        cores := 1,
        profile := rec(
        ),
        default_lang := ""
   )),

   ARMV7L := Class(ARMV7L, ARM, rec(
       cpuname := "ARMV7L",
       cores := 4,
   )),

   BlueGeneL := rec()
);


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

#
# SSE conversion bridges
#

# simple cases - signed<->unsigned

Class(_cvt_int_uint_saturated, ISA_Bridge_I, rec(
    props := ["saturation"],
    code  := (self, y, x, opts) >>
        assign( self._y(y,0), tcast(self.isa_to.t, bin_and(self._x(x,0), mask_gt(self._x(x,0), self.isa_from.t.zero()))) )
));

ISA_Bridge.add(Class(CVT_SSE_4x32ui_4x32i_wrap, ISA_Bridge_tcast, rec(
    isa_from := SSE_4x32f(T_Int(32)),   isa_to := SSE_4x32f(T_UInt(32)),
)));

ISA_Bridge.add(Class(CVT_SSE_4x32ui_4x32i_sat, _cvt_int_uint_saturated, rec(
    isa_from := SSE_4x32f(T_Int(32)),   isa_to := SSE_4x32f(T_UInt(32))
)));

ISA_Bridge.add(Class(CVT_SSE_4x32i_4x32ui_wrap, ISA_Bridge_tcast, rec(
    isa_from := SSE_4x32f(T_UInt(32)),  isa_to := SSE_4x32f(T_Int(32)),
)));

ISA_Bridge.add(Class(CVT_SSE_8x16ui_8x16i_wrap, ISA_Bridge_tcast, rec(
    isa_from := SSE_8x16i(T_Int(16)),   isa_to := SSE_8x16i(T_UInt(16)),
)));

ISA_Bridge.add(Class(CVT_SSE_8x16ui_8x16i_sat, _cvt_int_uint_saturated, rec(
    isa_from := SSE_8x16i(T_Int(16)),   isa_to := SSE_8x16i(T_UInt(16))
)));

ISA_Bridge.add(Class(CVT_SSE_8x16i_8x16ui_wrap, ISA_Bridge_tcast, rec(
    isa_from := SSE_8x16i(T_UInt(16)),  isa_to := SSE_8x16i(T_Int(16)),
)));

ISA_Bridge.add(Class(CVT_SSE_16x8ui_16x8i_wrap, ISA_Bridge_tcast, rec(
    isa_from := SSE_16x8i(T_Int(8)),    isa_to := SSE_16x8i(T_UInt(8)),
)));

ISA_Bridge.add(Class(CVT_SSE_16x8ui_16x8i_sat, _cvt_int_uint_saturated, rec(
    isa_from := SSE_16x8i(T_Int(8)),    isa_to := SSE_16x8i(T_UInt(8)),
)));

ISA_Bridge.add(Class(CVT_SSE_16x8i_16x8ui_wrap, ISA_Bridge_tcast, rec(
    isa_from := SSE_16x8i(T_UInt(8)),   isa_to := SSE_16x8i(T_Int(8)),
)));

# more complex cases

ISA_Bridge.add(Class(CVT_SSE_4x32i_4x32f_trunc, ISA_Bridge_I, rec(
    isa_from    := SSE_4x32f(T_Real(32)),
    isa_to      := SSE_4x32f(T_Int(32)),
    props       := ["trunc"],
    code := (self, y, x, opts) >> assign( self._y(y,0), vcvtt_4x32_f2i(self._x(x,0)) ),
)));

ISA_Bridge.add(Class(CVT_SSE_4x32i_4x32f_round, ISA_Bridge_I, rec(
    isa_from    := SSE_4x32f(T_Real(32)),
    isa_to      := SSE_4x32f(T_Int(32)),
    props       := ["round"],
    code := (self, y, x, opts) >> assign( self._y(y,0), vcvt_4x32_f2i(self._x(x,0)) ),
)));

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

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

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

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

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

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

ISA_Bridge.add(Class(CVT_SSE_8x16i_4x32i_sat, ISA_Bridge_I, rec(
    isa_from    := SSE_4x32f(T_Int(32)),
    isa_to      := SSE_8x16i(T_Int(16)),
    props       := ["saturation"],
    code := (self, y, x, opts) >> assign( self._y(y,0), vpacks_4x32i(self._x(x,0), self._x(x,1)) )
)));

# saturated SSE_4x32i to SSE_8x16ui
ISA_Bridge.add(Class(CVT_SSE_8x16ui_4x32i_sat, ISA_Bridge_I, rec(
    isa_from    := SSE_4x32f(T_Int(32)),
    isa_to      := SSE_8x16i(T_UInt(16)),
    props       := ["saturation"],
    code := (self, y, x, opts) >> let(
                c1 := self.isa_from.t.value(32768),
                c2 := self.isa_to.t.value(32768),
                u1 := var.fresh_t("U", self.isa_from.t),
                u2 := var.fresh_t("U", self.isa_from.t),
                m1 := var.fresh_t("U", self.isa_from.t),
                m2 := var.fresh_t("U", self.isa_from.t),
                decl([u1, u2, m1, m2], chain(
                    assign( u1, self._x(x,0) ),
                    assign( u2, self._x(x,1) ),
                    assign( m1, sub(bin_and(u1, mask_lt(u1.t.zero(), u1)), c1)),
                    assign( m2, sub(bin_and(u2, mask_lt(u1.t.zero(), u2)), c1)),
                    assign( self._y(y,0), add(tcast(self.isa_to.t, vpacks_4x32i(m1, m2)), c2) )
                ))
            ),
)));

# saturated SSE_4x32ui to SSE_8x16ui
ISA_Bridge.add(Class(CVT_SSE_8x16ui_4x32ui_sat, ISA_Bridge_I, rec(
    isa_from    := SSE_4x32f(T_UInt(32)),
    isa_to      := SSE_8x16i(T_UInt(16)),
    props       := ["saturation"],
    code := (self, y, x, opts) >> let(
                c1 := self.isa_from.t.value(32768),
                c2 := self.isa_to.t.value(32768),
                m1 := var.fresh_t("U", self.isa_from.t),
                m2 := var.fresh_t("U", self.isa_from.t),
                decl([m1, m2], chain(
                    assign( m1, sub(self._x(x,0), c1)),
                    assign( m2, sub(self._x(x,1), c1)),
                    assign( self._y(y,0), add(tcast(self.isa_to.t, vpacks_4x32i(m1, m2)), c2) )
                ))
            ),
)));

# SSE_4x32i to SSE_8x16ui without saturation by throwing out high word.
ISA_Bridge.add(Class(CVT_SSE_8x16ui_4x32i_wrap, ISA_Bridge_I, rec(
    isa_from    := SSE_4x32f(T_Int(32)),
    isa_to      := SSE_8x16i(T_UInt(16)),
    props       := ["wraparound"],
    code := (self, y, x, opts) >> let(
                u1 := var.fresh_t("U", self.isa_to.t),
                u2 := var.fresh_t("U", self.isa_to.t),
                m1 := var.fresh_t("U", self.isa_to.t),
                m2 := var.fresh_t("U", self.isa_to.t),
                m3 := var.fresh_t("U", self.isa_to.t),
                m4 := var.fresh_t("U", self.isa_to.t),
                decl([u1, u2, m1, m2, m3, m4], chain(
                    assign( u1, tcast(u1.t, self._x(x,0)) ),
                    assign( u2, tcast(u2.t, self._x(x,1)) ),
                    assign( m1, vunpacklo_8x16i(u1, u2) ),
                    assign( m2, vunpackhi_8x16i(u1, u2) ),
                    assign( m3, vunpacklo_8x16i(m1, m2) ),
                    assign( m4, vunpackhi_8x16i(m1, m2) ),
                    assign( self._y(y,0), vunpacklo_8x16i(m3, m4) )
                ))
            ),
)));

ISA_Bridge.add(Class(CVT_SSE_16x8i_8x16i_sat, ISA_Bridge_I, rec(
    isa_from    := SSE_8x16i(T_Int(16)),
    isa_to      := SSE_16x8i(T_Int(8)),
    props       := ["saturation"],
    code := (self, y, x, opts) >> assign( self._y(y,0), vpacks_8x16i(self._x(x,0), self._x(x,1)) )
)));

ISA_Bridge.add(Class(CVT_SSE_16x8ui_8x16i_sat, ISA_Bridge_I, rec(
    isa_from    := SSE_8x16i(T_Int(16)),
    isa_to      := SSE_16x8i(T_UInt(8)),
    props       := ["saturation"],
    code := (self, y, x, opts) >> assign( self._y(y,0), vpackus_8x16i(self._x(x,0), self._x(x,1)) )
)));

ISA_Bridge.add(Class(CVT_SSE_16x8ui_8x16ui_sat_as, ISA_Bridge_I, rec(
    isa_from    := SSE_8x16i(T_UInt(16)),
    isa_to      := SSE_16x8i(T_UInt(8)),
    props       := ["saturation"],
    code := (self, y, x, opts) >> let(
                t  := TVect(T_Int(16), 8),
                a1 := var.fresh_t("U", t),
                a2 := var.fresh_t("U", t),
                m1 := var.fresh_t("U", self.isa_to.t),
                m2 := var.fresh_t("U", self.isa_to.t),
                decl([a1,a2,m1,m2], chain(
                    assign( a1, arith_shr(tcast(t, self._x(x,0)), 15) ),
                    assign( a2, arith_shr(tcast(t, self._x(x,1)), 15) ),
                    assign( m1, vpacks_8x16i(a1, a2) ),
                    assign( m2, vpackus_8x16i(self._x(x,0), self._x(x,1)) ),
                    assign( self._y(y,0), bin_or(m1,m2) )
                ))
            ),
)));

ISA_Bridge.add(Class(CVT_SSE_16x8ui_8x16ui_sat_gt, ISA_Bridge_I, rec(
    isa_from    := SSE_8x16i(T_UInt(16)),
    isa_to      := SSE_16x8i(T_UInt(8)),
    props       := ["saturation"],
    code := (self, y, x, opts) >> let(
                t  := TVect(T_Int(16), 8),
                a1 := var.fresh_t("U", t),
                a2 := var.fresh_t("U", t),
                m1 := var.fresh_t("U", self.isa_to.t),
                m2 := var.fresh_t("U", self.isa_to.t),
                decl([a1,a2,m1,m2], chain(
                    assign( a1, gt(t.zero(), tcast(t, self._x(x,0))) ),
                    assign( a2, gt(t.zero(), tcast(t, self._x(x,1))) ),
                    assign( m1, vpacks_8x16i(a1, a2) ),
                    assign( m2, vpackus_8x16i(self._x(x,0), self._x(x,1)) ),
                    assign( self._y(y,0), bin_or(m1,m2) )
                ))
            ),
)));

ISA_Bridge.add(Class(CVT_SSE_8x16i_16x8i_gt, ISA_Bridge_I, rec(
    isa_from    := SSE_16x8i(T_Int(8)),
    isa_to      := SSE_8x16i(T_Int(16)),
    props       := [],
    code := (self, y, x, opts) >> let(
                sign := var.fresh_t("U", self.isa_from.t),
                decl( [sign], chain(
                    assign(sign, gt(self.isa_from.t.zero(), self._x(x,0))),
                    assign(self._y(y,0), tcast( self.isa_to.t, vunpacklo_16x8i(self._x(x,0), sign)) ),
                    assign(self._y(y,1), tcast( self.isa_to.t, vunpackhi_16x8i(self._x(x,0), sign)) )
                ))
            ),
)));

ISA_Bridge.add(Class(CVT_SSE_8x16i_16x8i_as, ISA_Bridge_I, rec(
    isa_from    := SSE_16x8i(T_Int(8)),
    isa_to      := SSE_8x16i(T_Int(16)),
    props       := [],
    code := (self, y, x, opts) >> chain(
                    assign(self._y(y,0), arith_shr(tcast( self.isa_to.t, vunpacklo_16x8i(self._x(x,0), self._x(x,0))), 8) ),
                    assign(self._y(y,1), arith_shr(tcast( self.isa_to.t, vunpackhi_16x8i(self._x(x,0), self._x(x,0))), 8) )
                ),
)));

ISA_Bridge.add(Class(CVT_SSE_8x16i_16x8ui, ISA_Bridge_I, rec(
    isa_from    := SSE_16x8i(T_UInt(8)),
    isa_to      := SSE_8x16i(T_Int(16)),
    props       := [],
    code := (self, y, x, opts) >> chain(
                    assign(self._y(y,0), vcastizxlo(self._x(x,0)) ),
                    assign(self._y(y,1), vcastizxhi(self._x(x,0)) )
                )
)));

ISA_Bridge.add(Class(CVT_SSE_4x32i_8x16i_gt, ISA_Bridge_I, rec(
    isa_from    := SSE_8x16i(T_Int(16)),
    isa_to      := SSE_4x32f(T_Int(32)),
    props       := [],
    code := (self, y, x, opts) >> let(
                sign := var.fresh_t("U", self.isa_from.t),
                decl( [sign], chain(
                    assign(sign, gt(self.isa_from.t.zero(), self._x(x,0))),
                    assign(self._y(y,0), tcast( self.isa_to.t, vunpacklo_8x16i(self._x(x,0), sign)) ),
                    assign(self._y(y,1), tcast( self.isa_to.t, vunpackhi_8x16i(self._x(x,0), sign)) )
                ))
            ),
)));

ISA_Bridge.add(Class(CVT_SSE_4x32i_8x16i_as, ISA_Bridge_I, rec(
    isa_from    := SSE_8x16i(T_Int(16)),
    isa_to      := SSE_4x32f(T_Int(32)),
    props       := [],
    code := (self, y, x, opts) >> chain(
                    assign(self._y(y,0), arith_shr(tcast( self.isa_to.t, vunpacklo_8x16i(self._x(x,0), self._x(x,0))), 16) ),
                    assign(self._y(y,1), arith_shr(tcast( self.isa_to.t, vunpackhi_8x16i(self._x(x,0), self._x(x,0))), 16) )
                )
)));

ISA_Bridge.add(Class(CVT_SSE_4x32i_8x16ui, ISA_Bridge_I, rec(
    isa_from    := SSE_8x16i(T_UInt(16)),
    isa_to      := SSE_4x32f(T_Int(32)),
    props       := [],
    code := (self, y, x, opts) >> chain(
                    assign(self._y(y,0), tcast( self.isa_to.t, vunpacklo_8x16i(self._x(x,0), self.isa_from.t.zero())) ),
                    assign(self._y(y,1), tcast( self.isa_to.t, vunpackhi_8x16i(self._x(x,0), self.isa_from.t.zero())) )
                )
)));

ISA_Bridge.add(Class(CVT_SSE_4x32f_4x32i, ISA_Bridge_I, rec(
    isa_from    := SSE_4x32f(T_Int(32)),
    isa_to      := SSE_4x32f(T_Real(32)),
    props       := [],
    code := (self, y, x, opts) >> assign(self._y(y,0), vcvt_4x32_i2f(self._x(x,0)))
)));