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

###############################################################################
##
#R  IsSelfSimRep
##
##  This is how IsSelfSim object is stored in GAP:
##  IsSelfSim object is a thing of kind "w = (w_1, w_2, ..., w_d)\pi", where
##    deg = d - arity of tree;
##    perm = \pi - permutation on first level;
##    w, w_1, ..., w_d - elements of free group representing elements of
##      automata group;
##    word = w;
##    states = [w_1, ..., w_d].
##
DeclareRepresentation("IsSelfSimRep",
                      IsComponentObjectRep and IsAttributeStoringRep,
                      ["word", "states", "perm", "deg"]);

InstallGlobalFunction(__AG_CreateSelfSim,
function(family, word, states, perm, invertible)
  local a, cat;

  if invertible then
    cat := IsInvertibleSelfSim and IsSelfSimRep;

    if perm^-1=fail then
      Error(perm, " is not invertible");
    else
      perm := AG_PermFromTransformation(perm);
    fi;
  else
    cat := IsSelfSim and IsSelfSimRep;
  fi;

  a := Objectify(NewType(family, cat),
                 rec(word := word,
                     states := states,
                     perm := perm,
                     deg := family!.deg));

  SetIsActingOnBinaryTree(a, a!.deg = 2);

  return a;
end);

###############################################################################
##
#M  SelfSim( <word>, <fam> )
##
InstallMethod(SelfSim, "for [IsAssocWord, IsSelfSimFamily]",
              [IsAssocWord, IsSelfSimFamily],
function(w, fam)
  local exp, wstates, curstate, newstate, curletter, newletter,
        nperm, i, j, perm, a, wtmp, reduced, invertible;

  if fam!.use_rws then
    w := AG_ReducedForm(fam!.rws, w);
  fi;

  if Length(w) = 0 then
    return One(fam);
  elif Length(w) = 1 then
    if ExponentSyllable(w, 1) = 1 then
      return fam!.recurgens[GeneratorSyllable(w, 1)];
    else
      if IsBound(fam!.recurgens[GeneratorSyllable(w, 1) + fam!.numstates]) then
        return fam!.recurgens[GeneratorSyllable(w, 1) + fam!.numstates];
      else
        Error(fam!.recurgens[GeneratorSyllable(w, 1)], " is not invertible");
      fi;
    fi;
  fi;

  exp := LetterRepAssocWord(w);
  for i in [1..Length(exp)] do
    if exp[i] < 0 then
      if IsBound(fam!.recurlist[-exp[i] + fam!.numstates]) then
        exp[i] := -exp[i] + fam!.numstates;
      else
        Error(fam!.recurgens[-exp[i]], " is not invertible");
      fi;
    fi;
  od;
  wstates := [];
  nperm := ();
  for i in [1..Length(exp)] do
    nperm := nperm * fam!.recurlist[exp[i]][fam!.deg+1];
  od;

  for i in [1..fam!.deg] do
    wstates[i] := One(fam!.freegroup);
    perm := ();

    for j in [1..Length(exp)] do
      wstates[i] := wstates[i]*fam!.recurgens[exp[j]]!.states[i^perm];
      perm := perm * fam!.recurlist[exp[j]][fam!.deg+1];
    od;

#    if fam!.use_rws and Length(wstates[i]) > 0 then
#      wstates[i] := AG_ReducedForm(fam!.rws, wstates[i]);
#    fi;
  od;

  invertible := true;
  if not fam!.isgroup then
    for i in exp do
      if i <= fam!.numstates and not IsInvertibleSelfSim(fam!.recurgens[i]) then
        invertible := false;
        break;
      fi;
    od;
  fi;

  return __AG_CreateSelfSim(fam, w, wstates, nperm, invertible);
end);

###############################################################################
##
#M  SelfSim( <word>, <a> )
##
InstallMethod(SelfSim, "for [IsAssocWord, IsSelfSim]", [IsAssocWord, IsSelfSim],
function(w, a)
  return SelfSim(w, FamilyObj(a));
end);

InstallMethod(MappedWord, [IsAssocWord,
                           IsList and IsAssocWordCollection,
                           IsList and IsSelfSimCollection],
function(w, fgens, agens)
  local img;
  img := MappedWord(w, fgens, List(agens, a -> a!.word));
  return SelfSim(img, FamilyObj(agens[1]));
end);

###############################################################################
##
#M  SelfSim( <word>, <list> )
##
InstallMethod(SelfSim, "for [IsAssocWord, IsList]",
                   [IsAssocWord, IsList],
function(w, list)
  local fam;
  fam := SelfSimFamily(list);
  if fam = fail then
    return fail;
  fi;
  return SelfSim(w, fam);
end);

###############################################################################
##
#M  PrintObj( <a> )
##
InstallMethod(PrintObj, "for [IsSelfSim]",
              [IsSelfSim],
function (a)
  local deg, printword, i;

  printword := function(w)
    if IsOne(w) then Print(AG_Globals.identity_symbol);
    else Print(w); fi;
  end;

  if true then
    View(a);
    return;
  fi;

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

###############################################################################
##
#M  ViewObj( <a> )
##
InstallMethod(ViewObj, "for [IsSelfSim]",
              [IsSelfSim],
function (a)
  if IsOne(a!.word) then Print(AG_Globals.identity_symbol);
  else Print(a!.word); fi;
end);

###############################################################################
##
#M  String( <a> )
##
InstallMethod(String, "for [IsSelfSim]",
              [IsSelfSim],
function (a)
  if IsOne(a!.word) then return AG_Globals.identity_symbol;
  else return String(a!.word); fi;
end);

###############################################################################
##
#M  Perm( <a> )
##
InstallMethod(Perm, "for [IsSelfSim]", [IsSelfSim],
function(a)
    return a!.perm;
end);

###############################################################################
##
#M  Word(<a>)
##
InstallMethod(Word, "for [IsSelfSim]", [IsSelfSim],
function(a)
    return a!.word;
end);

###############################################################################
##
#M  <a1> * <a2>
##
InstallMethod(\*, "for [IsSelfSim, IsSelfSim]", [IsSelfSim, IsSelfSim],
function(a1, a2)
    local a, i, fam, word, states;

    fam := FamilyObj(a1);
    word := a1!.word * a2!.word;

    if fam!.use_rws then
      word := AG_ReducedForm(fam!.rws, word);
    fi;

    if IsOne(word) then
      return One(a1);
    fi;

    states := List([1..a1!.deg], i -> a1!.states[i] * a2!.states[i^(a1!.perm)]);

    if fam!.use_rws then
      for i in [1..a1!.deg] do
        states[i] := AG_ReducedForm(fam!.rws, states[i]);
      od;
    fi;

    return __AG_CreateSelfSim(FamilyObj(a1), word, states, a1!.perm * a2!.perm,
                           IsInvertibleSelfSim(a1) and IsInvertibleSelfSim(a2));
end);

###############################################################################
##
#M  a1 < a2
##
InstallMethod(\<, "for [IsSelfSim, IsSelfSim]", IsIdenticalObj, [IsSelfSim, IsSelfSim],
function(a1, a2)
  local d, checked, checked_words, p, pos, np, np_words, i, perm1, perm2, cmp, G1, G2;

  G1 := GroupOfSelfSimFamily(FamilyObj(a1));
  G2 := GroupOfSelfSimFamily(FamilyObj(a2));

  # if a1 or a2 are elements of semigroup, which is not a group, then G1 or G2 is `fail'

  if IsIdenticalObj(G1, G2) and HasIsFinite(G1) and IsFinite(G1) and a1 = a2 then
    return false;
  fi;

  d := a1!.deg;

  checked := [[a1, a2]];
  checked_words := [[a1!.word, a2!.word]];
  pos := 0;

  while Length(checked) <> pos do
    pos := pos + 1;
    p := checked[pos];
    perm1 := p[1]!.perm;
    perm2 := p[2]!.perm;
    cmp := AG_TrCmp(perm1, perm2, d);
    if cmp < 0 then
      return true;
    elif cmp > 0 then
      return false;
    fi;
    for i in [1..d] do
      np := [Section(p[1], i), Section(p[2], i)];
      np_words := List(np, x -> x!.word);
      if not np_words in checked_words then
        Add(checked, np);
        Add(checked_words, np_words);
      fi;
    od;
  od;

  return false;
end);

###############################################################################
##
#M  InverseOp( <a> )
##
InstallMethod(InverseOp, "for [IsInvertibleSelfSim]", [IsInvertibleSelfSim],
function(a)
  local i, inv, fam, word, states;

  fam := FamilyObj(a);
  word := a!.word ^ -1;
  if fam!.use_rws then
    word := AG_ReducedForm(fam!.rws, word);
    if IsOne(word) then
      return One(a);
    fi;
  fi;

  states := List([1..a!.deg], i -> a!.states[i^(a!.perm^-1)]^-1);

  if fam!.use_rws then
    for i in [1..a!.deg] do
      states[i] := AG_ReducedForm(fam!.rws, states[i]);
    od;
  fi;

  return __AG_CreateSelfSim(FamilyObj(a), word, states, a!.perm^-1, true);
end);

###############################################################################
##
#M  OneOp( <a> )
##
InstallMethod(OneOp, "for [IsSelfSim]", [IsSelfSim],
function(a)
    return One(FamilyObj(a));
end);

###############################################################################
##
#M  StatesWords( <a> )
##
InstallMethod(StatesWords, "for [IsSelfSim]", [IsSelfSim],
function(a)
  return a!.states;
end);

###############################################################################
##
#M  Sections( <a> )
##
InstallMethod(Sections, "for [IsSelfSim]", [IsSelfSim],
function(a)
  return List(a!.states, s -> SelfSim(s, a));
end);

###############################################################################
##
#M  Section( <a>, <k>)
##
InstallMethod(Section, "for [IsSelfSim, IsPosInt]", [IsSelfSim, IsPosInt],
function(a, k)
  if k > a!.deg then
    Error("in Section(IsSelfSim, IsPosInt): invalid vertex ", k);
  fi;
  return SelfSim(a!.states[k], a);
end);

###############################################################################
##
#M  Section(a, seq)
##
## TODO
InstallMethod(Section, "for [IsSelfSim, IsList]", [IsSelfSim, IsList],
function(a, v)
  if Length(v) = 0 then
    return a;
  fi;

  if Length(v) = 1 then
    return Section(a, v[1]);
  fi;

  return Section(Section(a, v[1]), v{[2..Length(v)]});
end);

###############################################################################
##
#M  k ^ a
##
InstallMethod(\^, "for [IsPosInt, IsSelfSim]", [IsPosInt, IsSelfSim],
function(k, a)
    return k ^ Perm(a);
end);

###############################################################################
##
#M  seq ^ a
##
InstallMethod(\^, "for [IsList, IsSelfSim]", [IsList, IsSelfSim],
function(seq, a)
    local i, deg, img, cur;

    deg := DegreeOfTree(a);
    for i in seq do
      if not IsInt(i) or i < 1 or i > deg then
        Print("\^(IsList, IsFGSelfSim): ",
              i, "is out of range 1..", deg, "\n");
        return seq;
      fi;
    od;

    if Length(seq) = 0 then return []; fi;
    if Length(seq) = 1 then return [seq[1]^Perm(a)]; fi;

    return Concatenation([seq[1]^Perm(a)], seq{[2..Length(seq)]}^Section(a, seq[1]));
end);

###############################################################################
##
#M  PermOnLevelOp( <a>, <k>)
##
## TODO
InstallMethod(PermOnLevelOp, "for [IsInvertibleSelfSim, IsPosInt]",
              [IsInvertibleSelfSim, IsPosInt],
function(a, k)
  local dom, perm;

  if k = 1 then
    return a!.perm;
  fi;

  dom := AsList(Tuples([1..a!.deg], k));
  perm := List(dom, s -> s^a);
  perm := PermListList(dom, perm);

  return perm;
end);

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

###############################################################################
##
#M  IsActingOnBinaryTree( <a> )
##
InstallMethod(IsActingOnBinaryTree, "for [IsSelfSim]",
              [IsSelfSim],
function(a)
    return a!.deg = 2;
end);

###############################################################################
##
#M  IsOne( <a> )
##
InstallMethod(IsOne, "for [IsSelfSim]",
              [IsSelfSim],
function(a)
  local st, state_words, IsOneLocal, G;

  G := GroupOfSelfSimFamily(FamilyObj(a));

  if HasIsFinite(G) and IsFinite(G) then
    return IsOne(Image(IsomorphismPermGroup(G), a));
  fi;

  if HasIsContracting(G) and IsContracting(G) then
    return IsOne(Image(MonomorphismToAutomatonGroup(G), a));
  fi;

  IsOneLocal := function(s)
    local i, triv;
    if not IsOne(Perm(s)) then return false; fi;
    if s!.word in state_words then  return true; fi;

    Add(state_words, s!.word);
    triv := true;
    i := 1;
    while triv and i<=s!.deg do
      triv := IsOneLocal(Section(s, i));
      i := i+1;
    od;
    return triv;
  end;

  state_words := [];

  return IsOneLocal(a);
end);

###############################################################################
##
#M  a1 = a2
##
##
InstallMethod(\=, "for [IsSelfSim, IsSelfSim]", IsIdenticalObj, [IsSelfSim, IsSelfSim],
function(a1, a2)
  local areequalstates, d, checked_pairs, G1, G2;

  G1 := GroupOfSelfSimFamily(FamilyObj(a1));
  G2 := GroupOfSelfSimFamily(FamilyObj(a2));

  # if a1 or a2 are elements of semigroup, which is not a group, then G1 or G2 is `fail'

  if IsIdenticalObj(G1, G2) and HasIsFinite(G1) and IsFinite(G1) then
    return IsOne(Image(IsomorphismPermGroup(G1), a1*a2^-1));
  fi;

  d := a1!.deg;

  areequalstates := function(a, b)
    local j;

    if a!.word = b!.word then
      return true;
    fi;

    if [a!.word, b!.word] in checked_pairs then
      return true;
    else
      if a!.perm <> b!.perm then
        return false;
      fi;
      AddSet(checked_pairs, [a!.word, b!.word]);
      for j in [1..d] do
        if not areequalstates(Section(a, j), Section(b, j)) then
          return false;
        fi;
      od;
      return true;
    fi;
  end;

  checked_pairs := [];
  return areequalstates(a1, a2);
end);

InstallMethod(OrderUsingSections, "[IsSelfSim, IsCyclotomic]", true,
              [IsSelfSim, IsCyclotomic],
function(a, max_depth)
  local OrderUsingSections_LOCAL, cur_list, F, degs, vertex, AreConjugateUsingSmallRels, gens_ord2, CyclicallyReduce, res;

  CyclicallyReduce := function(w)
    local i, j, wtmp, reduced;

    for i in [1..Length(w)] do
      if -w[i] in gens_ord2 then w[i] := -w[i]; fi;
    od;

    repeat
      reduced := true;
      j := 1;
      while reduced  and j<Length(w) do
        if w[j]=-w[j+1] or (w[j]=w[j+1] and w[j] in gens_ord2) then
          reduced := false;
          wtmp := ShallowCopy(w{[1..j-1]});
          Append(wtmp, w{[j+2..Length(w)]});
          w := wtmp;
        fi;
        j := j+1;
      od;
    until reduced;

    repeat
      if Length(w)<2 then return w; fi;
      reduced := true;
      if w[1]=-w[Length(w)] or (w[1]=w[Length(w)] and w[1] in gens_ord2) then
        w := w{[2..Length(w)-1]};
        reduced := false;
      fi;
    until reduced;

    return w;
  end;

  AreConjugateUsingSmallRels := function(g, h)
    local i, g_list, h_list, long_cycle, l;
    g_list := CyclicallyReduce(LetterRepAssocWord(g));
    h_list := CyclicallyReduce(LetterRepAssocWord(h));
    if Length(g_list)<>Length(h_list) then return false; fi;
    l := [2..Length(g_list)];
    Add(l, 1);
    long_cycle := PermList(l);
    for i in [0..Length(g_list)-1] do
      if h_list=Permuted(g_list, long_cycle^i) then return true; fi;
    od;
    return false;
  end;

  OrderUsingSections_LOCAL := function(g)
    local i, el, orb, Orbs, res, st, reduced_word, loc_order;
    if IsOne(g) then return 1; fi;

    if IsActingOnBinaryTree(g) and IsSphericallyTransitive(g) then
      Info(InfoAutomGrp, 3, g!.word, " acts transitively on levels and is obtained from (", a!.word, ")^", Product(degs{[1..Length(degs)]}), "\n    by taking sections and cyclic reductions at vertex ", vertex);
      return infinity;
    fi;

    for i in [1..Length(cur_list)] do
      el := cur_list[i];
      if (AreConjugateUsingSmallRels(g!.word, el!.word) or AreConjugateUsingSmallRels((g!.word)^(-1), el!.word)) then
        if Product(degs{[i..Length(degs)]})>1 then
          if i>1 then Info(InfoAutomGrp, 3, "(", a!.word, ")^", Product(degs{[1..i-1]}), " has ", el!.word, " as a section at vertex ", vertex{[1..i-1]}); fi;
          Info(InfoAutomGrp, 3, "(", el!.word, ")^", Product(degs{[i..Length(degs)]}), " has conjugate of ", g!.word, " as a section at vertex ", vertex{[i..Length(degs)]});
          SetIsFinite(GroupOfSelfSimFamily(FamilyObj(a)), false);
          return infinity;
        else
#          Info(InfoAutomGrp, 3, "The group <G> might not be contracting, ", g, " has itself as a section.");
          return 1;
        fi;
      fi;
    od;
    if Length(cur_list)>=max_depth then return fail; fi;
    Add(cur_list, g);
    Orbs := OrbitsPerms([g!.perm], [1..g!.deg]);
    loc_order := 1;

    for orb in Orbs do
      Add(degs, Length(orb));
      Add(vertex, orb[1]);
      st := Section(g^Length(orb), orb[1]);
      reduced_word := AssocWordByLetterRep(FamilyObj(st!.word), CyclicallyReduce(LetterRepAssocWord(st!.word)));
#      Print(st!.word, " at ", vertex, "\n");
      res := OrderUsingSections_LOCAL(SelfSim(reduced_word, FamilyObj(g)));
      if res = infinity then return res;
      elif res=fail then
        loc_order:=fail;
      fi;
      if loc_order<>fail then
        loc_order := Lcm(loc_order, res*Length(orb));
      fi;
      Remove(degs);
      Remove(vertex);
    od;
    Remove(cur_list);
    return loc_order;
  end;

  F := FamilyObj(a)!.freegroup;
  if IsObviouslyFiniteState(FamilyObj(a)) then
    gens_ord2 := GeneratorsOfOrderTwo(FamilyObj(a));
  else
    gens_ord2 := [];
  fi;
  cur_list := [];
# degs traces at what positions we raise to what power
  degs := []; vertex := [];
  res := OrderUsingSections_LOCAL(a);
  if res=infinity then
    SetIsFinite(GroupOfSelfSimFamily(FamilyObj(a)), false);
    SetOrder(a, infinity);
  fi;
  return res;
end);

InstallMethod(OrderUsingSections, "for [IsSelfSim]", true,
              [IsSelfSim],
function(a)
  return OrderUsingSections(a, infinity);
end);

InstallMethod(SphericalIndex, "for [IsSelfSim]", [IsSelfSim],
function(a)
  # XXX check uses of SphericalIndex everywhere
  return rec(start := [], period := [a!.deg]);
end);

# XXX check uses of this everywhere
InstallMethod(DegreeOfTree, "for [IsSelfSim]", [IsSelfSim],
function(a)
  return a!.deg;
end);

# XXX check uses of this everywhere
InstallMethod(TopDegreeOfTree, "for [IsSelfSim]", [IsSelfSim],
function(a)
  return a!.deg;
end);

###############################################################################
##
#M  CanEasilyTestSphericalTransitivity( <a> )
##
InstallTrueMethod(CanEasilyTestSphericalTransitivity,
                  IsActingOnBinaryTree and IsSelfSim and IsFiniteState);

###############################################################################
##
#M  IsSphericallyTransitive( <a> )
##
InstallMethod(IsSphericallyTransitive, "for [IsInvertibleSelfSim]",
              [IsInvertibleSelfSim],
function(a)
  local w, i, ab, abs;

  if IsOne(Word(a)) then
    Info(InfoAutomGrp, 3, "IsSphericallyTransitive(a): false");
    Info(InfoAutomGrp, 3, "  IsOne(Word(a)): a = ", a);
    return false;
  fi;

  TryNextMethod();
end);

#########################################################################
##
#M  Order( <a> )
##
InstallMethod(Order, "for [IsInvertibleSelfSim]", true,
                   [IsInvertibleSelfSim],
function(a)
  local ord_loc;
  if HasIsFiniteState(GroupOfSelfSimFamily(FamilyObj(a))) and IsFiniteState(GroupOfSelfSimFamily(FamilyObj(a))) then
    if IsGeneratedByBoundedAutomaton(UnderlyingAutomatonGroup(GroupOfSelfSimFamily(FamilyObj(a)))) then
      return OrderUsingSections(a, infinity);
    fi;
  fi;
  if IsActingOnBinaryTree(a) and IsSphericallyTransitive(a) then
    return infinity;
  fi;
  ord_loc := OrderUsingSections(a, 10);
  if ord_loc <> fail then
    return ord_loc;
  fi;
  return OrderUsingSections(a, infinity);
end);

#########################################################################
##
#M  IsTransitiveOnLevel( <a>, <lev> )
##
InstallMethod(IsTransitiveOnLevel, "for [IsInvertibleSelfSim, IsPosInt]",
              [IsInvertibleSelfSim, IsPosInt],
function(a, lev)
  return Length(OrbitPerms([PermOnLevel(a, lev)], 1)) = a!.deg^lev;
end);

#########################################################################
##
#M  IsFiniteState( <a> )
##
InstallMethod(IsFiniteState, "for [IsSelfSim]", [IsSelfSim],
function(a)
  local st, state_words, find_all_sections;

  find_all_sections := function(s)
    local i;
    if not s!.word in state_words then
      Add(state_words, s!.word);
      for i in [1..s!.deg] do find_all_sections(Section(s, i)); od;
    fi;
  end;

  state_words := [];
  find_all_sections(a);
  SetAllSections(a, List(state_words, x -> SelfSim(x, FamilyObj(a))));
  return true;
end);

#########################################################################
##
#M  AllSections( <a> )
##
InstallMethod(AllSections, "for [IsSelfSim]", [IsSelfSim],
function(a)
  if IsFiniteState(a) then return AllSections(a); fi;
end);

#E

max_stars_repo_name: gap-packages/automgrp

max_issues_repo_name: gap-packages/automgrp

max_forks_repo_name: gap-packages/automgrp


text:
Length("abc");
# or same result with
Size("abc");

ext: gap

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
--/**
-- * <copyright>
-- *
-- * Copyright (c) 2008, 2009 IBM Corporation and others.
-- * All rights reserved.   This program and the accompanying materials
-- * are made available under the terms of the Eclipse Public License v2.0
-- * which accompanies this distribution, and is available at
-- * http://www.eclipse.org/legal/epl-v20.html
-- *
-- * Contributors:
-- *   IBM - Initial API and implementation
-- *   E.D.Willink - Lexer and Parser refactoring to support extensibility and flexible error handling
-- *   Adolfo Sanchez-Barbudo Herrera (Open Canarias) - LPG v 2.0.17 adoption (242153)
-- *   Adolfo Sanchez-Barbudo Herrera (Open Canarias) - Introducing new LPG templates (299396)
-- *
-- * </copyright>
-- */
--
-- The OCL Backtracking KeyWord Lexer, which is nominally identical to the
-- normal KeyWord Lexer, however the extra ERROR_TOKEN symbol makes it difficult
-- to share reliably.
--

%options slr
%options fp=OCLBacktrackingKWLexer,prefix=Char_
%options noserialize
%options package=org.eclipse.ocl.parser.backtracking
%options template=../../lpg/KeywordTemplateF.gi
%options export_terminals=("OCLBacktrackingParsersym.java", "TK_")
%options include_directory="..;../../lpg"

%Import
	OCLKWLexer.gi
%End

ext: gi


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

Class(HCOLCodegen, DefaultCodegen, rec(
    PointWise := (self, o, y, x, opts) >> let(i := Ind(o.N),
        loop(i, o.N, assign(nth(y, i), o.op.at(nth(x, i), i)))),

    ISumReduction := (self, o, y, x, opts) >> let(op := o.op, idval := o.idval, c:=o.child(1),
        n := Rows(c), t := TempVec(TArray(y.t.base_t(), n)), i := Ind(n), j := Ind(n),
        chain(
            loop(j, n,
                assign(nth(y,j), idval)),

            Cond(o.isSaturated = False,
                loop(o.var, o.domain, decl([t],
                    chain(
                    self(o.child(1), t, x, opts),
                    loop(i, n,
                          assign(nth(y,i), op(nth(y,i), nth(t,i))))
                ))),
                loopc(o.var, o.domain, ApplyFunc(logic_and, List([0..n-1], e->logic_neg(o.isSaturated(nth(y,e))))), decl([t],
                    chain(
                    self(o.child(1), t, x, opts),
                    loop(i, n,
                          assign(nth(y,i), op(nth(y,i), nth(t,i))))
                )))
            )
        )),

	COND := (self, o, y, x, opts) >> let(
			#Print("Cond:",o.cond, "\nTrue: ", o.child(1), "\nFalse:", o.child(2),"\n"),
			t := TempVec(TArray(TBool,1)),
			c1:= o.child(1), c2 := o.child(2),

			decl([t],
				chain(
					self(o.cond, t, x, opts),
					IF( nth(t, 0),
						self(o.child(1), y, x, opts),
						self(o.child(2), y, x, opts))))),

    Constant := (self, o, y, x, opts) >> assign(y, o.value),

    eT := (self, o, y, x, opts) >> assign(nth(y,0), nth(x, o.base)),

    eUnion := (self, o, y, x, opts) >> When(o.N=0 and o.base=-1, chain(), assign(nth(y, o.base), nth(x,0)) ),

    Mul1 := (self, o, y, x, opts) >> assign(nth(y,0), o.element*nth(x, 0)),

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

    ISumUnion := (self, o, y, x, opts) >>
       loop(o.var, o.domain, self(o.child(1), y, x, opts)),

    Inductor := (self, o, y, x, opts) >> assign(nth(y,0), cond(eq(o.i,0), o.initval, o.op.at(nth(y,0), nth(x,0)))),

    BinOp := (self, o, y, x, opts) >> let(i := Ind(o.N),
        loop(i, V(o.N), assign(nth(y,i), o.op.at(nth(x, i), nth(x,i+o.N))))),

    OLCompose := meth(self, o,y,x,opts)
        local ch, numch, vecs;
        ch    := o.children();
        numch := Length(ch);

        # propagate x and y arrays and create temporary arrays
        vecs  := self._composePropagate(ch, y, x, c -> TempArray(y,x,c));

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

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

));

max_stars_repo_name: spiral-software/spiral-package-hcol

max_issues_repo_name: spiral-software/spiral-package-hcol

max_forks_repo_name: spiral-software/spiral-package-hcol


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

# Auxiliary Functions for Symbols and Non-Terminals
# =================================================
# MP, from 08/17/00

#F Stride Permutations
#F -------------------
#F

# Let K a field and d | n. Then the stride permutation L_d^n is a
# permutation of degree n defined by
#   L_d^n (e_i tensor e_j) = (e_j tensor e_i)
#   for all i, j
# where e_i is the ith canonical base vector in K^d
# and e_j the jth in K^(n/d).
# Hence L_d^n exchanges the tensor components of dimension
# d and n/d resp.
# L_d^n can also be viewed the following way. Let x be
# a (d x n/d)-matrix stored in row major order. L_d^n permutes
# x such that it is stored in column major order
# (i.e. transposes x).

# Lemma: L_d^n is given as a permutation of 0..n-1 by
#   i -> i * d mod (n - 1), 0 <= i < n-1.
#
# Lemma: The following computation rules hold:
#   a) L_r^rst * L_s^rst = L_rs^rst.
#   b) L_t^rst = (L_t^rt tensor 1_s)(1_r tensor L_t^st).
#
# Lemma: L_2^(2i) has exactly 2 fixpoints.
#
# Lemma: (has to be proven)
#   p prime: ord(L_(p^i)^(p^n)) = n/gcd(n, i).

#F StridePerm( <n>, <d> )
#F   returns the stride permutation L_d^n (see above)
#F

StridePerm := function ( n, d )
  local M;

  # check arguments
  if not ( IsInt(n) and IsInt(d) and n > 0 and d > 0 ) then
    Error("<n> and <d> must be positive integers");
  fi;
  if not n mod d = 0 then
    Error("<d> must divide <n>");
  fi;

  # create matrix and transpose
  M := List([1..d], i -> [(i-1)*n/d+1..i*n/d]);

  return
    MappingPermListList(
      Concatenation(TransposedMat(M)),
      [1..n]
    );
end;

# Computes the same thing, but slightly slower.

StridePerm1 := function ( n, d )
  local L;

  # check arguments
  if not ( IsInt(n) and IsInt(d) and n > 0 and d > 0 ) then
    Error("<n> and <d> must be positive integers");
  fi;
  if not n mod d = 0 then
    Error("<d> must divide <n>");
  fi;

  # compute images according to lemma above
  L    := List([1..n], i -> (i - 1) * d mod (n - 1) + 1);
  L[n] := n;

  # return permutation
  return MappingPermListList([1..n], L);
end;

# CycleTypeStridePerm( <n>, <d> )
#   computes the cycle type of L_d^n. The cycle type is a
#   list of pairs [ .., [c_i, n_i], .. ], which denotes
#   that the permutation contains exactly n_i cycles of
#   length c_i.
#   If c_i = 1, then n_i is the number of fixpoints.
#

CycleTypeStridePerm := function ( n, d )
  return Collected(CycleLengths(StridePerm(n, d), [1..n]));
end;

#F Twiddle Matrices
#F ----------------
#F

# The twiddle factors appear in the Cooley-Tukey FFT, namely,
# for d|n we have:
#   F_n = (F_n/d tensor I_d) * T(n, d) * (I_n/d tensor F_d) * L(n, n/d)
# where F_n is a DFT of size n, and L denotes a stride permutation
# defined as above

#F TwiddleDiag( <n>, <d> [, <k> ] )
#F TwiddleDiag( <list-of-n-d[-k]> )
#F   returns the twiddle factors for <d>|<n> as a list:
#F     T(n, d) = direct_sum_(i = 1)^n/d diag(w_n^0, ..., w_n^(d-1)^i
#F   If the optional parameter <k> is given then w_n is replaced by
#F   w_n^k, gcd(<n>, <k>) = 1.
#F

TwiddleDiag := function ( arg )
  local n, d, k;

  if Length(arg) = 1 and IsList(arg[1]) then
    arg := arg[1];
  fi;

  if Length(arg) = 2 then
    n := arg[1];
    d := arg[2];
    k := 1;
  elif Length(arg) = 3 then
    n := arg[1];
    d := arg[2];
    k := arg[3];
  else
    Error(
      "usage:\n",
      "  TwiddleDiag( <n>, <d> [, <k> ] )\n",
      "  TwiddleDiag( <list-of-n-d[-k]> )"
    );
  fi;

  # check arguments
  if not ( IsInt(n) and IsInt(d) and n > 0 and d > 0 ) then
    Error("<n> and <d> must be positive integers");
  fi;
  if not n mod d = 0 then
    Error("<d> must divide <n>");
  fi;
  if not Gcd(n, k) = 1 then
    Error("gcd(<n>, <k>) must be 1");
  fi;

  return
    Concatenation(
      List([1..n/d], i -> List([1..d], j -> E(n)^(k*(j-1)*(i-1))))
    );
end;


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

Class(DTTBase, NonTerminal, rec(
    abbrevs := [ N -> Checked(IsInt(N), N >= 1, [N]) ],
    dims := self >> [self.params[1], self.params[1]],
    isReal := True,

    normalizedArithCost := self >> let(n:=Double(self.params[1]),
        2.5*n*d_log(n)/d_log(2)),

    hashAs := self >> self
));

Declare(DCT3, DCT7);

#F DCT1(<n>) - Discrete Cosine Transform, Type I, non-terminal
#F Definition: (n x n)-matrix [ cos(k*l*pi/(n-1)) | k,l = 0...n-1 ]
#F Note:       The natural size for a DCT1 is 2^k + 1
#F Example:    DCT1(9)
#F
Class(DCT1, DTTBase, rec(
    terminate := self >> Mat(DCT_Iunscaled(self.params[1])),
    transpose := self >> Copy(self),
    SmallRandom := () -> Random([3,4,5,7,9,13,17,25,33])
));

#F DCT2(<n>) - Discrete Cosine Transform, Type II, non-terminal
#F Definition: (n x n)-matrix [ cos(k*(l+1/2)*pi/n) | k,l = 0...n-1 ]
#F Note:       DCT2 is the transpose of DCT3
#F Example:    DCT2(8)
Class(DCT2, DTTBase, rec(
    terminate := self >> Mat(DCT_IIunscaled(self.params[1])),
    transpose := self >> DCT3(self.params[1]),
    SmallRandom := () -> Random([2,3,4,5,6,8,9,10,12,15,16,18,24,27,30,32])
));

#F DCT3(<n>) - Discrete Cosine Transform, Type III, non-terminal (unscaled)
#F Definition: (n x n)-matrix [ cos((k+1/2)*l*pi/n) | k,l = 0...n-1 ]
#F  [scaled]   (n x n)-matrix [ a_l*cos((k+1/2)*l*pi/n) | k,l = 0...n-1 ]
#F                              a_j = 1/sqrt(2) for j = 0 and = 1 else
#F Note:       DCT3 is the transpose of DCT2, scaled NOT supported yet
#F Example:    DCT3(8)
Class(DCT3, DTTBase, rec(
    terminate := self >> Mat(DCT_IIIunscaled(self.params[1])),
    transpose := self >> DCT2(self.params[1]),
    SmallRandom := () -> Random([2,3,4,5,6,8,9,10,12,15,16,18,24,27,30,32])
));

#F DCT4(<size>) - Discrete Cosine Transform, Type IV, non-terminal
#F Definition: (n x n)-matrix  [ cos((k+1/2)*(l+1/2)*pi/n) | k,l = 0...n-1 ]
#F Note:       DCT4 is symmetric
#F Example:    DCT4(8)
Class(DCT4, DTTBase, rec(
    terminate := self >> Mat(DCT_IVunscaled(self.params[1])),
    transpose := self >> Copy(self),
    SmallRandom := () -> Random([2, 3, 4, 6, 8, 9, 12, 16, 18, 24, 27, 32])
));

#F DCT5(<size>) - Discrete Cosine Transform, Type V, non-terminal
#F Definition: (n x n)-matrix   [ cos(k*l*pi/(n-1/2)) | k,l = 0...n-1 ]
#F Note:       DCT5 is symmetric
#F Example:    DCT5(8)
Class(DCT5, DTTBase, rec(
    terminate := self >> Mat(DCT_Vunscaled(self.params[1])),
    transpose := self >> Copy(self),
));

#F DCT6(<size>) - Discrete Cosine Transform, Type VI, non-terminal
#F Definition: (n x n)-matrix   [ cos(k*(l+1/2)*pi/(n-1/2)) | k,l = 0...n-1 ]
#F Note:       The transpose of DCT6 is DCT7
#F Example:    DCT6(8)
Class(DCT6, DTTBase, rec(
    terminate := self >> Mat(DCT_VIunscaled(self.params[1])),
    transpose := self >> DCT7(self.params[1]),
));

#F DCT7(<size>) - Discrete Cosine Transform, Type VII, non-terminal
#F Definition: (n x n)-matrix   [ cos((k+1/2)*l*pi/(n-1/2)) | k,l = 0...n-1 ]
#F Note:       The transpose of DCT7 is DCT6
#F Example:    DCT7(8)
Class(DCT7, DTTBase, rec(
    terminate := self >> Mat(DCT_VIIunscaled(self.params[1])),
    transpose := self >> DCT6(self.params[1]),
));

#F DCT8(<size>) - Discrete Cosine Transform, Type VIII, non-terminal
#F Definition: (n x n)-matrix [ cos((k+1/2)*(l+1/2)*pi/(n+1/2)) | k,l = 0...n-1 ]
#F Note:       DCT8 is symmetric
#F Example:    DCT8(8)
Class(DCT8, DTTBase, rec(
    terminate := self >> Mat(DCT_VIIIunscaled(self.params[1])),
    transpose := self >> Copy(self)
));


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

Class(VPRDFT, PRDFT);

PRF12_VCT_Children := (vlen, N,k,PRFt,DFTt,PRFtp,PRF1) -> let(m:=vlen, n:=N/vlen,
	[[ PRFt(m,k), DFTt(m,k), PRF1(n,k), PRFtp(m,k) ]]
);

Class(CompileBarrier, Buf);
CompileBarrier.code := (self, y, x) >> chain(skip(), compiler._CodeSums(self.child(1), y, x), skip());

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

    (SUM(
	RC(Scat(H(Nf+1,mf+1,0,n))) * C[1] * Gath(H(2*(nf+1)*m, m, 0, 1)), #2*(nf+1))),

	When(nc=1, [],
	ISum(j,
	     RC(Scat(BH(Nf+1,N,m,j+1,n))) *
	     Conj * RC(C[2]) * Tw(j) * L(2*m, m) *
	     RC(Gath(H((nf+1)*m, m, (j+1)*m, 1))))),  #j+1, nf+1

	When(IsOddInt(n), [],
	RC(Scat(H(Nf+1,mc,nf,n))) * C[4] * Gath(H(2*(nf+1)*m, m, 2*m*nf, 1)))#2*nf, 2*(nf+1))))
    )) *
    #L(Rows(C[3])*m,m)) *
    CompileBarrier((transforms.VTensor(C[3], m)))
);

RulesFor(VPRDFT, rec(
    #F PRDFT1_CT: projection of DFT_CT
    VPRDFT_VCT := rec(
	isApplicable := P -> not IsPrime(P[1]),
	vlen := 4,
	allChildren  := (self, P) >> PRF12_VCT_Children(self.vlen, P[1], P[2], PRDFT1, DFT1, PRDFT3, PRDFT1),
	rule := (P,C) -> let(N:=P[1], k:=P[2], m:=Cols(C[1]),
	    PRF12_VCT_Rule(N, k, C, Diag(BHD(m,1,-1)), j->RC(Diag(fPrecompute(Twid(N,m,k,0,0,j+1)))))))
));

# opts := CopyFields(PRDFTDefaults, rec(globalUnrolling:=0, unparser:=CMacroUnparser, includes:=["<include/real_sse.h>"], compflags:="-O1 -msse -fomit-frame-pointer"));
# c:=CodeRuleTreeOpts(RandomRuleTree(VPRDFT(16)), opts);
# PrintCode("sub1", c, opts);
# PropagateTypes(c);
# CMacroUnparser(c, 4, 4);

# r := RandomRuleTree(VPRDFT(64));
# s := RecurseOpts(r, PRDFTDefaults);
# s1 := Collect(s, RecursStep)[2];
Import(search);

VPRDFT_example_gcc := function(n)
    SwitchRulesByName(PRDFT, [PRDFT1_Base1, PRDFT1_Base2, PRDFT1_CT, PRDFT_Rader, PRDFT_PD, PRDFT1_PF]);
    SwitchRulesByName(PRDFT3, [PRDFT3_Base1, PRDFT3_Base2, PRDFT3_CT, PRDFT3_OddToPRDFT1]);

    Local.opts := CopyFields(spiral.PRDFTDefaults,
	rec(globalUnrolling:=34,
	    unparser:=compiler.CMacroUnparser,
	    includes:=["<include/real_sse.h>"], compflags:="-O3 -msse2 -fomit-frame-pointer -ffast-math"));
    if not IsBound(Local.hash) then Local.hash := HashTableDP(); fi;

    return DP(VPRDFT(n), rec(hashTable := Local.hash), Local.opts);
end;

VPRDFT_example := function(n)
    SwitchRulesByName(PRDFT, [PRDFT1_Base1, PRDFT1_Base2, PRDFT1_CT, PRDFT_Rader, PRDFT_PD, PRDFT1_PF]);
    SwitchRulesByName(PRDFT3, [PRDFT3_Base1, PRDFT3_Base2, PRDFT3_CT, PRDFT3_OddToPRDFT1]);

    Local.opts := CopyFields(spiral.PRDFTDefaults,
	rec(globalUnrolling:=34,
	    unparser:=compiler.CMacroUnparser,
	    includes:=["<include/real_sse.h>"], compflags:="-O3"));
    if not IsBound(Local.hash) then Local.hash := HashTableDP(); fi;

    return DP(VPRDFT(n), rec(hashTable := Local.hash), Local.opts);
end;


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

Class(PDTT_Base, TaggedNonTerminal, rec(
    abbrevs := [ n -> Checked(IsPosIntSym(n), [n]) ],
    dims := self >> [self.params[1], self.params[1]],
    isReal := True,
    print := (self,i,is) >> Print(self.__name__, "(", PrintCS(self.params), ")",
	When(self.transposed, Print(".transpose()")),
        When(self.tags<>[], Print(".withTags(", self.tags, ")"))),

    normalizedArithCost := self >> let(n := self.params[1],
	floor(2.5 * n * log(n) / log(2.0)) )
));

#
# Discrete Trigonometric Transforms (DTTs) via PRF
#
Declare(PDST3, PDCT3);

Class(PDST4, PDTT_Base, rec(
    transpose := self >> Copy(self),
    terminate := self >> Mat(DST_IVunscaled(EvalScalar(self.params[1]))),
));

Class(PDCT4, PDTT_Base, rec(
    transpose := self >> Copy(self),
    terminate := self >> Mat(DCT_IVunscaled(EvalScalar(self.params[1]))),
));

Class(PDCT2, PDTT_Base, rec(
    terminate := self >> Mat(DCT_IIunscaled(EvalScalar(self.params[1]))),
    transpose := self >> PDCT3(self.params[1])
));

Class(PDST2, PDTT_Base, rec(
    terminate := self >> Mat(DST_IIunscaled(EvalScalar(self.params[1]))),
    transpose := self >> PDST3(self.params[1])
));

Class(PDCT3, PDTT_Base, rec(
    terminate := self >> Mat(DCT_IIIunscaled(EvalScalar(self.params[1]))),
    transpose := self >> PDCT2(self.params[1])
));

Class(PDST3, PDTT_Base, rec(
    terminate := self >> Mat(DST_IIIunscaled(EvalScalar(self.params[1]))),
    transpose := self >> PDST2(self.params[1])
));

RulesFor(PDST4, rec(
    PDST4_Base2 := DST4_Base, #BaseRule(PDST4, 2),
    PDST4_CT := rec(
	isApplicable     := P -> P[1] > 2, #not IsPrime(P),
	forTransposition := false,
	allChildren  := P -> let(N := P[1], Map2(DivisorPairs(2*N), (m,n) -> Cond(
		IsEvenInt(m) and IsEvenInt(n), [ PRDFT3(m,-1).transpose(), PRDFT3(n) ],
		#IsEvenInt(m) and IsEvenInt(n), [ PDHT3(m).transpose(), PDHT3(n) ],
		IsEvenInt(m) and IsOddInt(n),  [ PRDFT3(m,-1).transpose(), PRDFT3(n), PDST4(m/2) ],
		IsOddInt(m)  and IsEvenInt(n), [ PRDFT3(m,-1).transpose(), PRDFT3(n), PDST4(n/2) ],
		IsOddInt(m)  and IsOddInt(n),  Error("This can't happen")))),

	rule := (P,C) -> let(N := P[1], m := Rows(C[1]), n := Cols(C[2]),
	    nf:=Int(n/2), nc:=Int((n+1)/2), mf:=Int(m/2), mc:=Int((m+1)/2), j:=Ind(nf), i:=Ind(mf),
	    # mult twid by -E(4) for RDFT, by 1/2 for DHT
	    t := ScaledTwid(2*N, mf, 1, 1/2, 1/2, j, -E(4)),
	    T := When(IsOddInt(m), Diag(diagDirsum(t,fConst(1,1))), Diag(t)),
	    Et := Tensor(I(mf),Mat([[0,-1],[-1,0]])),
	    Ett := Tensor(I(mf),Mat([[0,-1],[1,0]])),
	    SUM(
		ISum(j,
		    Scat(BH(N, 2*N-1, m, j, n)) * C[1] * #Ett *
		    RC(T) *
		    #Et * # DHT only
		    RC(Gath(H(mc*nc, mc, j, nc)))
		),
		When(IsEvenInt(n), [], Scat(H(N,m/2,nf,n))*C[3]*Gath(H(2*mc*nc, m/2, 2*nf, 2*nc)))
	    ) *
	    SUM(
		ISum(i, RC(Scat(H(mc*nc, nc, i*nc, 1))) * C[2] * Gath(BH(N, 2*N-1, n, i, m))),
                # output is imaginary, but we output into real slot, because subseq transform
		# is jIR2, which we implement as IR2 * (-j)
		When(IsEvenInt(m), [], Scat(H(2*mc*nc, nc, 2*mf*nc, 2))*C[3]*Gath(H(N,n/2,mf,m)))
	    ))),

    PDST4_CT_SPL := rec(
	isApplicable     := P -> P[1] > 2 and ForAny(DivisorPairs(2*P[1]), d->IsEvenInt(d[1]) and IsEvenInt(d[2])),
	forTransposition := false,

	allChildren  := P -> let(N := P[1],
	    Map2(Filtered(DivisorPairs(2*N), d -> IsEvenInt(d[1]) and IsEvenInt(d[2])),
		(m,n) -> [ PRDFT3(m,-1).transpose(), PRDFT3(n) ])),

	rule := (P,C) -> let(N := P[1], m := Rows(C[1]), n := Cols(C[2]),
            TT := Diag(fPrecompute(diagMul(fConst(N/2, -E(4)),
                        fCompose(dOmega(8 * N, 1),
                            diagTensor(dLin(N/m, 2, 1, TInt), dLin(m/2, 2, 1, TInt)))))),

	    Prm(Refl(N, 2*N-1, N, L(m*n, m))) *
	    Tensor(I(n/2), C[1]) *
            RC(TT) *
            RC(L(m*n/4, n/2)) *
	    Tensor(I(m/2), C[2]) *
	    Prm(Refl(N, 2*N-1, N, L(m*n, m)))
	)
    )
));

RulesFor(PDCT4, rec(
    PDCT4_Base2 := rec(isApplicable:=P->P[1]=2, rule:=(P,C)->DCT4(2).terminate()),
    PDCT4_Base4 := rec(isApplicable:=P->P[1]=4,
	               allChildren := P -> [[ DCT4(4) ]],
		       rule:=(P,C)->C[1]),

    PDCT4_CT := rec(
	isApplicable     := P -> P[1] > 2, #not IsPrime(P),
	forTransposition := false,
	allChildren  := P -> let(N := P[1], Map2(DivisorPairs(2*N), (m,n) -> Cond(
		IsEvenInt(m) and IsEvenInt(n), [ PRDFT3(m,-1).transpose(), PRDFT3(n) ],
		#IsEvenInt(m) and IsEvenInt(n), [ PDHT3(m).transpose(), PDHT3(n) ],
		IsEvenInt(m) and IsOddInt(n),  [ PRDFT3(m,-1).transpose(), PRDFT3(n), PDCT4(m/2) ],
		IsOddInt(m)  and IsEvenInt(n), [ PRDFT3(m,-1).transpose(), PRDFT3(n), PDCT4(n/2) ],
		IsOddInt(m)  and IsOddInt(n),  Error("This can't happen")))),

	rule := (P,C) -> let(N := P[1], m := Rows(C[1]), n := Cols(C[2]),
	    nf:=Int(n/2), nc:=Int((n+1)/2), mf:=Int(m/2), mc:=Int((m+1)/2), j:=Ind(nf), i:=Ind(mf),
	    t := fPrecompute(Twid(2*N, mf, 1, 1/2, 1/2, j)),
	    T := When(IsOddInt(m), Diag(diagDirsum(t,fConst(1,1))), Diag(t)),
	    Et := Tensor(I(mf),Mat([[0,1],[1,0]])),
	    SUM(
		ISum(j,
		    Scat(BH(N, 2*N-1, m, j, n)) * Diag(BHN(m)) * C[1] *
		    # mult twid by 1/2 for DHT
		    RC(T) *
		    #Et * # DHT only
		    RC(Gath(H(mc*nc, mc, j, nc)))
		),
		When(IsEvenInt(n), [], Scat(H(N,m/2,nf,n))*C[3]*Gath(H(2*mc*nc, m/2, 2*nf, 2*nc)))
	    ) *
	    SUM(
		ISum(i, RC(Scat(H(mc*nc, nc, i*nc, 1))) * C[2] * Diag(BHN(n)) * Gath(BH(N, 2*N-1, n, i, m))),
		When(IsEvenInt(m), [], Scat(H(2*mc*nc, nc, 2*mf*nc, 2))*C[3]*Gath(H(N,n/2,mf,m)))
	    ))),
));

NewRulesFor(PDCT4, rec(
    PDCT4_CT_SPL := rec(
        libApplicable := t -> eq(imod(t.params[1], 2), 0),
	applicable := t -> IsSymbolic(t.params[1]) or (t.params[1] > 2 and (t.params[1] mod 4) = 0),
        extraLeftTags := [],
	forTransposition := false,

        freedoms := t -> [ divisorsIntNonTriv(t.params[1]/2) ],

	child  := (self, t, fr) >> let(N := t.params[1], f := fr[1],
	    [ spiral.sym.ASP(-1).RDFT3(2*f).withTags(self.extraLeftTags).transpose(),
	      spiral.sym.ASP.RDFT3(div(N,f)) ]),

	apply := (t,C,Nonterms) -> let(N := t.params[1], m := Rows(C[1]), n := Cols(C[2]),
            mh := Rows(C[1])/2, nh := Cols(C[2])/2,
            D := RC(Diag(fPrecompute(fCompose(dOmega(8 * N, 1),
                            diagTensor(dLin(N/m, 2, 1, TInt), dLin(m/2, 2, 1, TInt)))))),

            Grp(Scat(Refl1(nh, m)) *
                Tensor(I(nh), Diag(BHN(m)) * C[1]) *
                D *
                Tensor(I(nh), Tr(2, mh))) *
            #RC(Tr(mh, nh)) *
            Grp(Tr(mh, n) *
                Tensor(I(mh), C[2] * Diag(BHN(n))) *
                Gath(Refl1(mh, n)))
	)
    )
));

RulesFor(PDCT4, rec(
    # Variant better suited for vectorization
    PDCT4_CT_SPL_Vec := rec(
	isApplicable     := P -> P[1] > 2 and (P[1] mod 2) = 0,
	forTransposition := false,

	allChildren  := P -> let(N := P[1],
	    Map2(Filtered(DivisorPairs(2*N), d -> IsEvenInt(d[1]) and IsEvenInt(d[2])),
		(m,n) -> [ PRDFT3(m,-1).transpose(), PRDFT3(n) ])),

	rule := (P,C) -> let(N := P[1], m := Rows(C[1]), n := Cols(C[2]),
            TT := Diag(fPrecompute(fCompose(dOmega(8 * N, 1),
                        diagTensor(dLin(N/m, 2, 1, TInt), dLin(m/2, 2, 1, TInt))))),

            Prm(condIJ(N, n/2)) *
            Prm(L(N, m)) *
	    Tensor(I(n/2),
                   condM(m,m/2) * Diag(BHN(m)) * C[1]) *
            RC(TT) *
            Prm(L(N, n/2)) *
	    Tensor(I(m/2),
                   L(n, 2) * C[2] * Diag(BHN(n)) * condK(n, 2)) *
            Prm(L(N, m/2)) *
            Prm(condIJ(N, m/2))
	)
    )

));

testPDCT4 := function(n)
    local opts, r, s;
    opts := CopyFields(SpiralDefaults, rec(
            breakdownRules := rec(
                PDCT4 := [PDCT4_Base2, PDCT4_CT_SPL_Vec],
                PRDFT1 := [PRDFT1_Base1, PRDFT1_Base2, PRDFT1_CT],
                PRDFT3 := [PRDFT3_Base1, PRDFT3_Base2, PRDFT3_CT],
                DFT := [DFT_Base, DFT_CT]
            )));
    r := RandomRuleTree(PDCT4(n), opts);
    s := SumsRuleTree(r, opts);
    return [opts,r,s];
end;

NewRulesFor(PDCT2, rec(
    PDCT2_Base2 := rec(
	forTransposition := true,
        applicable := t -> t.params[1] = 2,
        apply := (t, C, Nonterms) -> DCT2(2).terminate()
    ),

    PDCT2_Base4 := rec(
	forTransposition := true,
        applicable := t -> t.params[1] = 4,
        apply := (t, C, Nonterms) ->
            LIJ(4) *
            DirectSum(
                Diag(FList(TReal, [ 1, 0.70710678118654757 ])) * F(2),
                Rot(cospi(13/8), sinpi(13/8)) * J(2)) *
            (Tensor(I(2), F(2))) ^ LIJ(4)
    ),

    PDCT2_CT_SPL := rec(
	forTransposition := true,
        libApplicable := t -> eq(imod(t.params[1], 2), 0),
        extraLeftTags := [],

	applicable := t -> let(N:=t.params[1],
            IsSymbolic(N) or
            (N > 2 and IsEvenInt(N) and (N mod 4) = 0)),

	freedoms := t -> [ divisorsIntNonTriv(t.params[1]/2) ],

        # N/2 = k*m
        child := (self, t, fr) >> let(N:=t.params[1], k:=fr[1], m:=N/2/k,
            [ PDCT2(2*k).withTags(self.extraLeftTags),
              spiral.sym.ASP(-1).RDFT3(2*k).transpose().withTags(self.extraLeftTags),
              spiral.sym.ASP(+1).URDFT(2*m) ]),

	apply := (t,C,Nonterms) -> let(N:=t.params[1], k:=Rows(C[1])/2, m:=Cols(C[3])/2,
            D := RC(Diag(fPrecompute(fCompose(dOmega(4 * N, 1),
                        diagTensor(dLin(m-1, 1, 1, TInt), dLin(k, 2, 1, TInt)))))),
            Grp(
                Scat(Refl0_u(m, 2*k)) *
                DirectSum(
                    C[1] * KK(k, 2),
                    Tensor(I(m-1), Diag(BHN(2*k)) * C[2]) * D
                )) *
            Grp(
                RC(Tr(k, m)) *
                Tensor(I(k), C[3]) *
                Gath(Refl1(k, 2*m))
            )
	)
    )
));

NewRulesFor(PDST2, rec(
    PDST2_Base2 := rec(
        applicable := t -> t.params[1] = 2,
        apply := (t, C, Nonterms) -> DST2(2).terminate()
    ),

    PDST2_CT_SPL := rec(
	forTransposition := true,
	applicable := t -> let(N:=t.params[1], N > 2 and IsEvenInt(N) and (N mod 4) = 0),
	freedoms := t -> [ DivisorsIntNonTriv(t.params[1]/2) ],

        # N/2 = k*m
        child := (t, fr) -> let(N:=t.params, k:=fr[1], m:=N/2/k,
                [ DST2(2*k), PRDFT3(2*k,-1).transpose(), URDFT(2*m)  ]),

	apply := (t,C,Nonterms) -> let(N:=t.params[1], k:=Rows(C[1])/2, m:=Cols(C[3])/2,
            TT := Diag(fPrecompute(fCompose(dOmega(4 * N, 1),
                        # ie multiply all twiddles by -E(4)
                        diagAdd(diagDirsum(fConst(k, 0), fConst(k*(m-1), 3*N)),
                                diagTensor(dLin(m, 1, 0, TInt), dLin(k, 2, 1, TInt)))))),

            J(N) *
            Kp(N, 2*k) *
	    DirectSum(
                J(2*k) * C[1] * Diag(BHN(2*k)) * K(2*k, 2),
                Tensor(I(m-1), J(2*k) * M(2*k, k) * C[2])
            ) *
            RC(TT) *
            RC(L(N/2, m)) *
	    Tensor(I(k), C[3] * Diag(BHN(2*m))) *
	    Prm(Refl(N, 2*N-1, N, L(2*N, 2*k)))
	)
    )

));

RulesFor(DCT5, rec(
  DCT5_Rader := rec(
	forTransposition := false,
	isApplicable     := P -> P > 2 and IsPrime(2*P-1),
	allChildren      := P -> [[ IPRDFT(P-1,-1), PRDFT(P-1,-1) ]],

	diag := N -> Sublist(
	    DFT_Rader.raderDiag(2*N-1,1,PrimitiveRootMod(2*N-1)),
	    [1..Int((N-1)/2)]*2),

	# 3rd col with 0's is for PRDFT, not necessary for (non-packed) RDFT
	raderMid := (self, N) >> let(Fsize := 2*N-1,
	             DirectSum(Mat([[1, 1, 0], [1, -1/(Fsize-1), 0], [0,0,0]]),
			       RC(Diag(FData(self.diag(N)))))),

	# NOTE, special case for even P-1
	rule := (self,P,C) >> let(N := P,
	    #RealRR(N).transpose() *

	    Gath(H(2*P-1, P, 0, 1)) *
	    RR(2*P-1).transpose() *
	    DirectSum(I(1), VStack(I(P-1), I(P-1))) *

	    DirectSum(I(1), C[1]) *
#		C[1]*Diag(1,1, Replicate(2*Int((P-2)/2), 2),
#		               When(IsEvenInt(P-1), [1,1],[]))) *
	    self.raderMid(N) *
	    DirectSum(I(1), C[2]) *

            # replace by RealRR
	    Gath(H(2*P-1, P, 0, 1)) *
	    RR(2*P-1) * DirectSum(I(1), VStack(I(P-1), J(P-1)))
    ))
));

#q:=ExpandSPL(DCT5(6))[1];
#s:=SPLRuleTree(q);
#me := MatSPL(s);
#them := remat(MatSPL(q.node));

#RDFT-11,  DCT5(6) 15(6/9 rots) + 19(12/7 irdft) + 17(12/5 rdft) = 30/21 = 51
#           DST5(5)
#           10(blk)
# fftw 60/50

# RDFT-13
#  DCT5(7)  18(8/10=4/8+0/1+2/1 rots) + 18(14/4 irdft) + 18(14/4 rdft) = 36/18=54
#  DST5(6)
#  12(blk)
#  fftw 76/34=110

# DFT-13
# fftw 176+68=244
# us = 54*2 (dct5) + 24(blk)*2 + dst5


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

--
-- 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 /.keywordKind[%rule_number] = ./

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

    $BeginAction /.%Header./

    $EndAction /../

    $BeginJava /.%BeginAction./

    $EndJava /.%EndAction./
%End

%Globals
    /.

    ./
%End

%Headers
    /.
    export class %action_type extends %prs_type
    {
        private  inputBytes : byte[];
        private final number keywordKind[] = new number[%num_rules + 1];

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

        public  lexer(curtok : number, lasttok : number) : number
        {
            number current_kind = getKind(inputBytes[curtok]),
                act;

            for (act = tAction(START_STATE, current_kind);
                 act > NUM_RULES && act < ACCEPT_ACTION;
                 act = tAction(act, current_kind))
            {
                curtok++;
                current_kind = (curtok > lasttok
                                       ? %eof_char
                                       : getKind(inputBytes[curtok]));
            }

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

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

        public void setInputBytes(byte[] inputBytes) { this.inputBytes = inputBytes; }

    ./
%End

%Rules
    /.

        public %action_type(byte[] inputBytes, number identifierKind)
        {
            this.inputBytes = inputBytes;
            keywordKind[0] = identifierKind;
    ./
%End

%Trailers
    /.

            for (i : number = 0; i < keywordKind.length; i++)
            {
                if (keywordKind[i] == 0)
                    keywordKind[i] = identifierKind;
            }
        }
    }
    ./
%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

_computeNumDen := (value) -> let(
    list := Filtered([1..10000], e -> IntDouble(value * e) = value * e),

    When(value = 0.0, [0, 1],
	 When(Length(list) = 0, [0, 1],
	      [IntDouble(list[1] * value), list[1]]))
);

_shiftElements := (sizeIdx, size, shft) -> let(
    numDen := _computeNumDen(shft), # numDen[1]/numDen[2] = shft
    i := Ind([0..sizeIdx - 1]),
    onlyPositive := (sizeIdx < size),
    # Lambda(i, cond(lt(i, size / 2), E(numDen[2] * size) ^ (sign(2 * i) * numDen[1] * i), E(numDen[2] * size) ^ (sign(2 * i - size) * numDen[1] * (i - size))))
    # Lambda(i, cond(lt(i, size / 2), E(numDen[2] * size) ^ (numDen[1] * i), E(numDen[2] * size) ^ (numDen[1] * (i - size))))
    Lambda(i, cond(onlyPositive,
                   E(numDen[2] * size) ^ (numDen[1] * i),
                   cond(lt(i, size / 2),
                        E(numDen[2] * size) ^ (numDen[1] * i),
                        E(numDen[2] * size) ^ (numDen[1] * (i - size)))))
);

#_tabelize := f->let(vals := f.tolist(), i := Ind(f.domain()), When(
#    ForAll(vals, v->v=vals[1]), Lambda(i, vals[1]),
#    Lambda(i, FData(List(vals, j->Value(TComplex, _unwrap(j)))).at(i))
#));

_tabelize := f->let(vals := f.tolist(), i := Ind(f.domain()), When(
    ForAll(vals, v->v=vals[1]), Lambda(i, vals[1]),
    let(    rtab := FData(Flat(List(vals, e->[re(e), im(e)]))),
        Lambda(i, cxpack(rtab.at(2*i), rtab.at(2*i+1))))
));

NewRulesFor(TResample, rec(
		   TResample_TGath := rec(
		       forTransposition := false,
		       switch           := true,
		       applicable       := nt -> ForAny(Zip2(nt.params[1], nt.params[2]), l-> l[1]<l[2]),
		       children      := nt -> [[TResample(nt.params[1], nt.params[1], nt.params[3]), TGath(_toBox(nt.params[2], nt.params[1]))]],
		       apply := (nt, C, cnt) -> C[1] * C[2]
		   ),

		   TResample_TScat := rec(
		       forTransposition := false,
		       switch           := true,
		       applicable       := nt -> ForAny(Zip2(nt.params[1], nt.params[2]), l-> l[1]>l[2]),
		       children      := nt -> [[TScat(_toBox(nt.params[1], nt.params[2])), TResample(nt.params[2], nt.params[2], nt.params[3])]],
		       apply := (nt, C, cnt) -> C[1] * C[2]

		   ),

		   TResample_MD_nofrac := rec(
		       forTransposition := false,
		       switch           := true,
		       applicable       := nt -> (nt.params[1] = nt.params[2]) and (Length(nt.params[1]) > 1) and ForAll(nt.params[3], e -> e = 0.0),
		       children := nt -> [[I(Product(nt.params[1]))]],
		       apply := (nt, C, cnt) -> C[1]
		   ),

		   # ASSUMPTION IS THAT FASTEST DIMENSION IS REAL DFT-ED
		   TResample_MD_frac := rec(
		       forTransposition := false,
		       switch           := true,
		       applicable       := nt -> (nt.params[1] = nt.params[2]) and (Length(nt.params[1]) > 1) and ForAny(nt.params[3], e -> (abs(e) > 0.0) and (abs(e) < 1.0)),
		       children      := nt -> [[IMDPRDFT(nt.params[1], 1), MDPRDFT(nt.params[2], -1)]],

		       apply := (nt, C, cnt) -> let(
    			   sizeX := Last(nt.params[2]),
    			   sizeXf := sizeX / 2 + 1,
    			   shftX := Last(nt.params[3]),

    			   dXf := diagMul(fConst(sizeXf, 1/sizeX), _shiftElements(sizeXf, sizeX, shftX)),

    			   sizes := DropLast(nt.params[2], 1),
    			   shfts := DropLast(nt.params[3], 1),

    			   df := List([1..Length(sizes)], e -> diagMul(fConst(sizes[e], 1/sizes[e]), _shiftElements(sizes[e], sizes[e], shfts[e]))),
                   dd := List(df::[dXf], _tabelize),

    			   C[1] *
    			   RCDiag(RCData(diagTensor(dd))) *
    			   C[2]
		       )
		   )
	       )
	   );

max_stars_repo_name: franzfranchetti/spiral-package-fftx

max_issues_repo_name: franzfranchetti/spiral-package-fftx

max_forks_repo_name: franzfranchetti/spiral-package-fftx


text:
# The Loop function will apply some function to every tuple built by taking
# the i-th element of each list. If one of them is exhausted before the others,
# the loop continues at its begining. Only the longests lists will be precessed only once.
Loop := function(a, f)
    local i, j, m, n, v;
    n := Length(a);
    v := List(a, Length);
    m := Maximum(v);
    for j in [1 .. m] do
        f(List([1 .. n], i -> a[i][1 + RemInt(j - 1, v[i])]));
    od;
end;

# Here we simply print each "row"
f := function(u)
    Perform(u, Print);
    Print("\n");
end;

Loop([["a", "b", "c"], ["A", "B", "C"], [1, 2, 3]], f);

aA1
bB2
cC3

Loop([["a", "b"], ["A", "B", "C", "D", "E"], [1, 2, 3]], f);

aA1
bB2
aC3
bD1
aE2

max_stars_repo_name: djgoku/RosettaCodeData

max_issues_repo_name: djgoku/RosettaCodeData

max_forks_repo_name: djgoku/RosettaCodeData


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

#F ==========================================================================
#F AParSMP(<num_threads>, <tid>) - SMP parallelization tag
#F    .params[1] becomes num_threads
#F    .params[2] becomes tid
#F
Class(AParSMP, AGenericTag, rec(
    isSMP := true,
    updateParams := meth(self)
        if Length(self.params)=1 then self.params := [self.params[1], threadId()];
        elif Length(self.params)=2 then ;
        else Error("Usage: AParSMP(<num_threads>, [<tid>])");
        fi;
    end
));

#F ==========================================================================
#F SMPBarrier(<nthreads>, <tid>, <spl>)
#F
#F ==========================================================================
Class(SMPBarrier, Buf, BaseContainer, rec(
    # inheriting from Buf helps to get extra methods from Buf, for example
    # .cannotChangeDataFormat (from paradigms/vector). Buf is considered a general wrapper.
    doNotMarkBB := true,
    new := (self, nthreads, tid, spl) >> Checked(IsPosInt0Sym(tid), IsPosIntSym(nthreads), IsSPL(spl),
        SPL(WithBases(self, rec(_children:=[spl], tid := tid, nthreads := nthreads, dimensions := spl.dims())))),

    dims := self >> self._children[1].dims(),

    rChildren := self >> [self.nthreads, self.tid, self._children[1]],
    rSetChild := meth(self, n, newC)
        if n=1 then self.nthreads := newC;
        elif n=2 then self.tid := newC;
        elif n=3 then self._children[1] := newC;
        else Error("<n> must be in [1..3]");
        fi;
    end
));

#F ==========================================================================
#F SMPSum(<nthreads>, <tid>, <var>, <domain>, <spl>) - parallel loop with <nthreads> threads
#F
#F As a matrix SMPSum(p, var, domain, spl) == ISum(var, domain, spl)
#F ==========================================================================
Class(SMPSum, ISum, rec(
    doNotMarkBB := true,
    abbrevs := [ (p, var, domain, spl) -> Checked(
                    IsInt(p) or IsScalar(p), IsVar(var), IsInt(p) or IsScalar(domain), IsSPL(spl),
                    [p, threadId(), var, domain, spl]) ],

    new := meth(self, nthreads, tid, var, domain, spl)
        local res;
        Constraint(IsSPL(spl));
	# if domain is an integer (not symbolic) it must be non-zero
	Constraint(IsPosIntSym(domain));
        var.isLoopIndex := true;
        var.range := domain;
	res := SPL(WithBases(self, rec(nthreads:=nthreads, tid:=tid, _children := [spl], var := var, domain := domain)));
	res.dimensions := res.dims();
	return res;
    end,

    rChildren := self >> [self.nthreads, self.tid, self.var, self.domain, self._children[1]],

    from_rChildren := (self, rch) >> ApplyFunc(ObjId(self), rch),

    rSetChild := meth(self, n, newChild)
	if n=1 then self.nthreads := newChild;
	elif n=2 then self.tid := newChild;
	elif n=3 then self.var := newChild;
	elif n=4 then self.var.range := newChild; self.var.isLoopIndex := true;
	elif n=5 then self._children[1] := newChild;
        else Error("<n> must be in [1..5]");
	fi;
    end,

    print := (self, i, is) >> Print(self.__name__, "(", self.nthreads, ", ", self.tid, ", ",
        self.var, ", ", self.domain, ",\n", Blanks(i+is), self.child(1).print(i+is, is), "\n",
	Blanks(i), ")", self.printA())
));

#F ==========================================================================
#F SMP(<num_threads>, <tid>, <spl>) - parallel container for rewriting
#F
#F SMP(ISum(...)) -- becomes --> SMPSum(...)
#F ==========================================================================
Class(SMP, BaseContainer, rec(
    doNotMarkBB := true,

    abbrevs := [ (nthreads, spl) -> Checked(IsInt(nthreads) or IsScalar(nthreads), IsSPL(spl),
                                            [nthreads, threadId(), spl]) ],

    new := (self, nthreads, tid, spl) >>
	SPL(WithBases(self, rec(nthreads:=nthreads, tid:=tid, dimensions := spl.dimensions, _children := [spl]))),

    rChildren := self >> [self.nthreads, self.tid, self._children[1]],

    rSetChild := meth(self, n, newChild)
        if   n = 1 then self.nthreads := newChild;
        elif n = 2 then self.tid := newChild;
	elif n = 3 then self._children[1] := newChild;
	else Error("<n> must be in [1..3]"); fi;
	return self;
    end,

    sums := self >> ObjId(self)(self.nthreads, self.tid, self.child(1).sums()),

    normalizedArithCost := self >> self.child(1).normalizedArithCost()
));


text:
roots := n -> List([0 .. n-1], k -> E(n)^k);

r:=roots(7);
# [ 1, E(7), E(7)^2, E(7)^3, E(7)^4, E(7)^5, E(7)^6 ]

List(r, x -> x^7);
# [ 1, 1, 1, 1, 1, 1, 1 ]


text:
MitM_CookieCount := 0;
MitM_CookieJar := rec();

InstallGlobalFunction(MitM_CompleteProcedure,
function(result, attr...)
    if Length(attr) = 0 then
        attr := rec();
    elif Length(attr) = 1 then
        attr := attr[1];
    else
        ErrorNoReturn("MitM_CompleteProcedure: takes 1 or 2 arguments (not ",
                      Length(attr) + 1, ")");
    fi;
    return MitM_OMRecToXML(MitM_OMRecToOMOBJRec(
                 OMATTR( attr
                       , OMA( OMS( "scscp1"
                                 , "procedure_completed" )
                            , result ) ) ) );
end);

InstallGlobalFunction(MitM_TerminateProcedure,
function(message, attr...)
    if Length(attr) = 0 then
        attr := rec();
    elif Length(attr) = 1 then
        attr := attr[1];
    else
        ErrorNoReturn("MitM_TerminateProcedure: takes 1 or 2 arguments (not ",
                      Length(attr) + 1, ")");
    fi;
    return MitM_OMRecToXML(MitM_OMRecToOMOBJRec(
                 OMATTR( attr
                       , OMA( OMS( "scscp1"
                                 , "procedure_terminated" )
                            , OME( OMS( "scscp1"
                                      , "error_system_specific" )
                                 , [ OMSTR(message) ] ) ) ) ) );
end);

InstallValue(MitM_SCSCPHandlers, rec(
    procedure_call := function(attr, oma)
        local t, eval, rattr, res;

        Info(InfoMitMServer, 15, " Evaluating... ", oma);
        Info(InfoMitMServer, 15, " Attributes... ", attr);

        t := NanosecondsSinceEpoch();
        eval := MitM_OMRecToGAP(oma);
        t := NanosecondsSinceEpoch() - t;

        if eval.success then
            rattr := rec( call_id := attr.call_id
                        , info_runtime := t / 1000000. );
            if IsBound(attr.option_return_cookie) then
                MitM_CookieJar.(MitM_CookieCount) := eval.result;
                res := MitM_GAPToOMRec(MitM_CookieCount);
                MitM_CookieCount := MitM_CookieCount + 1;
            elif IsBound(attr.option_return_object) then
                res := MitM_GAPToOMRec(eval.result);
            elif IsBound(attr.option_return_nothing) then
                res := rec();
            fi;

            return MitM_CompleteProcedure(res, rattr);
        else
            Info(InfoMitMServer, 15, " Error during evaluation: ", eval.error);
            return MitM_TerminateProcedure(eval.error, rec(call_id := attr.call_id));
        fi;
    end
) );

InstallGlobalFunction(MitM_HandleSCSCP,
function(node)
    local attr, scscp_call, scscp_oma, content;

    content := MitM_Content(node);
    # Validate wrt SCSCP v1.3 spec - procedure call (4.1.1)
    if not (Length(content) = 1 and MitM_Tag(content[1]) = "OMATTR") then
        Info(InfoMitMServer, 15, " Invalid procedure call: OMATTR expected");
        return MitM_TerminateProcedure("procedure call: OMOBJ should contain one OMATTR and nothing else");
    elif MitM_Tag(MitM_Content(content[1])[2]) <> "OMA" then
        Info(InfoMitMServer, 15, " Invalid procedure call: OMA expected");
        return MitM_TerminateProcedure("procedure call: OMOBJ: OMATTR's 2nd object should be an OMA");
    elif not (Length(MitM_Content(MitM_Content(content[1])[2])) = 2 and
              MitM_Tag(MitM_Content(MitM_Content(content[1])[2])[1]) = "OMS") then
        Info(InfoMitMServer, 15, " Invalid procedure call: OMS expected");
        return MitM_TerminateProcedure("procedure call: OMOBJ: OMATTR: OMA's 1st object should be an OMS");
    elif MitM_Tag(MitM_Content(MitM_Content(content[1])[2])[2]) <> "OMA" then
        Info(InfoMitMServer, 15, " Invalid procedure call: OMA expected");
        return MitM_TerminateProcedure("procedure call: OMOBJ: OMATTR: OMA's 2nd object should be an OMA");
    fi;

    attr := MitM_ATPToRec(MitM_Content(content[1])[1]);
    scscp_call := MitM_Content(MitM_Content(content[1])[2])[1];
    scscp_oma := MitM_Content(MitM_Content(content[1])[2])[2];

    if MitM_CD(scscp_call) = "scscp1" then
        if MitM_Name(scscp_call) = "procedure_call" then
            return MitM_SCSCPHandlers.procedure_call(attr, scscp_oma);
        elif MitM_Name(scscp_call) = "terminate_procedure" then
            # TODO:
        fi;
    else
        Info(InfoMitMServer, 15, " Unsupported CD ", scscp_call.attributes.cd);
        return fail;
    fi;
end);

InstallGlobalFunction(MitM_SCSCPHandler,
function(addr, stream)
    local done, version, r, reply, obj;
    done := false;

    Info(InfoMitMServer, 5, "Accepted MitM Connection on ", TCP_AddrToString(addr));
    version := MitM_SCSCPServerHandshake(stream, stream);
    Info(InfoMitMServer, 5, " SCSCP Protocol Version ", version);
    while not done do
        r := MitM_ReadSCSCP(stream);
        Info(InfoMitMServer, 15, " Received: ", r);
        if r.success <> true then
            Info(InfoMitMServer, 15, " Bad object received");
            Info(InfoMitMServer, 15, "  error: ", r.error);
            Info(InfoMitMServer, 15, " closing connection.");
            CloseStream(stream);
            done := true;
            # WriteLine(stream, Concatenation("error: ", r.error));
        else
            reply := MitM_HandleSCSCP(r.result);
            Info(InfoMitMServer, 15, " Evaluated to ", reply);
            WriteLine(stream, "<?scscp start ?>");
            WriteLine(stream, reply);
            WriteLine(stream, "<?scscp end ?>");
        fi;
    od;
    Info(InfoMitMServer, 5, "Leaving handler for ", TCP_AddrToString(addr));
end);

InstallGlobalFunction(StartMitMServer,
function(args...)
    local opt;

    opt := ShallowCopy(MitM_DefaultServerOptions);
    if IsBound(args[1]) and IsString(args[1]) then
        opt.hostname := args[1];
    fi;
    if IsBound(args[2]) and IsPosInt(args[2]) then
        opt.port := args[2];
    fi;

    Info(InfoMitMServer, 5, "Starting MitM TCP Server on ", opt.hostname, ":", opt.port);
    StartTCPServer(opt.hostname, opt.port, MitM_SCSCPHandler);
end);

InstallGlobalFunction(StreamToMitMServer,
function(hostname, port...)
    local stream;
    if Length(port) = 0 then
        port := 26133; # SCSCP default
    elif Length(port) = 1 then
        port := port[1];
    else
        Error("MitM_ConnectionToServer: 1 or 2 arguments expected, but ",
              Length(port) + 1, " found");
    fi;
    stream := ConnectInputOutputTCPStream(hostname, port);
    MitM_SCSCPClientHandshake(stream, stream);
    return stream;
end);

InstallGlobalFunction(SendObjToMitMServer,
function(stream, obj)
    local xml;
    xml := MitM_OMRecToXML(obj);
    WriteLine(stream, "<?scscp start ?>");
    WriteLine(stream, xml);
    WriteLine(stream, "<?scscp end ?>");
    return MitM_ReadSCSCP(stream);
end);

InstallGlobalFunction(MitM_ProcedureCall,
function(obj)
    local base64, call_id, opts, call, attr;
    # obj should be an OMA to be run
    base64 := Concatenation("ABCDEFGHIJKLMNOPQRSTUVWXYZ",
                            "abcdefghijklmnopqrstuvwxyz",
                            "0123456789+/");
    call_id := Concatenation(MitM_cdbase, ":",
                             List([1..16], i -> Random(base64)));
    # TODO: support option_return_cookie and option_return_nothing
    opts := rec(call_id := OMSTR(call_id), option_return_object := OMSTR(""));
    call := OMA(OMS("scscp1", "procedure_call"), obj);
    attr := OMATTR(opts, call);
    return OMOBJ(attr);
end);

InstallGlobalFunction(GetAllowedHeads,
function(stream)
    local call, out, obj, attr, oma_outer, oma_inner, list;
    call := MitM_ProcedureCall(OMA(OMS("scscp2", "get_allowed_heads")));
    out := SendObjToMitMServer(stream, call);
    if not out.success then
        return fail;
    fi;
    obj := out.result;
    # TODO: check that the server successfully sent a list of heads
    attr := MitM_Content(obj)[1];
    oma_outer := MitM_Content(attr)[2];
    oma_inner := MitM_Content(oma_outer)[2];
    list := MitM_Content(oma_inner);
    return list{[2..Length(list)]};
end);

max_stars_repo_name: markuspf/MathsInTheMiddle

max_issues_repo_name: markuspf/MathsInTheMiddle

max_forks_repo_name: markuspf/MathsInTheMiddle