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

# ================================================================
# Sum of Monomials
# NOTE: fix sum_mons, mon_diags, mon_perms
# ================================================================
sum_mons := pat -> pat.target(SUMAcc).cond(e->ForAll(e.children(), IS(Mon)));
mon_vdiags := pat -> List(pat.val.children(), c -> c.vdiag);
mon_sdiags := pat -> List(pat.val.children(), c -> c.sdiag);
mon_perms := pat -> List(pat.val.children(), c -> c.perm);

compat_tensor_mons := (p,v,s) -> let(
    pp := Filtered(p, c->domain(c) > 1),
    vv := Filtered(v, c->domain(c) > 1),
    ss := Filtered(s, c->domain(c) > 1),
    fully_compat_tensor_chains(pp, vv, domain, domain) and
    fully_compat_tensor_chains(vv, ss, domain, domain));

merge_tensor_mons := function(pp, vv, ss, combine) #, fid, gid, hid)
   local i, j, k, np, nv, ns, p, v, s, res;
   np := Length(pp);
   nv := Length(vv);
   ns := Length(ss);
   res := []; i:=1; j:=1; k := 1;

   while i <= np or j <= nv or k <= ns do
       if (i <= np and j <= nv and k <= ns) and
           domain(pp[i]) = domain(vv[j]) and
           domain(vv[j]) = domain(ss[k])
       then
           Add(res, combine(pp[i], vv[j], ss[k]));
           i := i+1;
           j := j+1;
           k := k+1;
       else
           p := fId(1);
           v := II(1);
           s := II(1);

           if i<=np and domain(pp[i]) = 1 then
               p := pp[i]; i:=i+1; fi;
           if j<=nv and domain(vv[j]) = 1 then
               v := vv[j]; j:=j+1; fi;
           if k<=ns and domain(ss[k]) = 1 then
               s := ss[k]; k:=k+1; fi;

           Add(res, combine(p,v,s));
       fi;
   od;
   return res;
end;

#UseRuleSet(RulesSums);

# pulling in sum of monomials into ISum
ARule( Compose,
       [ @(1, ISum), sum_mons(@(2)) ],
  e -> [ ISum(@(1).val.var, @(1).val.domain, @(1).val.child(1) * @(2).val).attrs(@(1).val) ]);

# Gath * SUMAcc(Mon, Mon,...) -> SUMAcc(GMon, GMon, ...)
ARule( Compose,
       [ [Gath, @(1)], sum_mons(@(2)) ],
  e -> [ SUMAcc(
	  let(vdiags := mon_vdiags(@(2)),
	      sdiags := mon_sdiags(@(2)),
	      perms  := mon_perms(@(2)),
	      numdiags := Length(vdiags),

	      List([1..numdiags],
		  i -> GMon(fCompose(perms[i], @(1).val),
		            fCompose(vdiags[i], perms[i], @(1).val),
		            fCompose(sdiags[i], perms[i], @(1).val))))) ]
);

# low-priority suck in of Gath into SUMAcc
ARule( Compose,
       [ @(1, SUMAcc), @(2, Gath) ],
  e -> [ ApplyFunc(ObjId(@(1).val),
	           List(@(1).val.children(), c -> c * @(2).val)) ]);

# Combining GMon and Gath, SMon and Scat, Mon and Prm
ARule( Compose,
       [ @(1, GMon), @(2, Gath) ], # o 1-> 2->
  e -> [ GMon(fCompose(@(2).val.func, @(1).val.perm),
	      @(1).val.vdiag, @(1).val.sdiag) ] );

ARule( Compose,
       [ @(1, Scat), @(2, SMon) ], # o 1-> 2->
  e -> [ GMon(fCompose(@(1).val.func, @(2).val.perm),
	      @(2).val.vdiag, @(2).val.sdiag) ] );

ARule( Compose,
       [ @(1, Mon), @(2, Prm) ], # o 1-> 2->
  e -> [ GMon(fCompose(@(2).val.func, @(1).val.perm),
	      fCompose(@(1).val.vdiag, @(2).val.func.transpose()),
	      fCompose(@(1).val.sdiag, @(2).val.func.transpose())) ]);

#UseRuleSet(RulesFuncSimp);

# ================================================================
# Simplification of Intervals
# ================================================================

full_interval := (int, size) ->
    [int.target(II), size, zero, @.cond(e->e=size.val)];

left_interval := (int, size, endd) ->
    [int.target(II), size, zero, endd];

right_interval := (int, size,start) ->
    [int.target(II), size, start, @.cond(e->e=size.val)];

# full interval o index mapping func @(3) = interval with domain of @(3)
# II(n, 0, n) o f = II(domain(f), 0, domain(f))
ARule( fCompose,
      [ full_interval(@, @(1)), @(2) ],
 e -> [ II(pdomain(@(2))) ]
);

# partial interval o dirsum of perms
# II(n, 0, k) o fDirsum(perm_k, perm_?) -> II(n,0,k)
ARule( fCompose,
      [ left_interval(@(0), @(1), @(2)),
	[fDirsum, @(3).cond(e -> is_perm(e) and range(e)=@(2).val), ...] ],
 e -> [ @(0).val ]
);

# partial interval o dirsum of perms
# II(n, n-k, n) o fDirsum(perm_?, perm_k) -> II(n,n-k,n)
ARule( fCompose,
      [ right_interval(@(0), @(1), @(2)),
	[fDirsum, ..., @(4).cond(e -> is_perm(e) and range(e) = @(1).val-@(2).val)] ],
 e -> [ @(0).val ]
);

# fTensor(IIn, IIk) -> II(n*k)
ARule( fTensor, [full_interval(@(0), @(1)), full_interval(@(0), @(2))],
 e -> [ II(@(1).val * @(2).val) ]
);

# fDirsum(IIn, IIk) -> II(n+k)
ARule( fDirsum, [full_interval(@(0), @(1)), full_interval(@(0), @(2))],
 e -> [ II(@(1).val + @(2).val) ]
);

# Interval Shifting, we support shift left (shift <= start), and right (shift >= end)
# II o Z = shifted II#
ARule( fCompose,
      [ [II, @(1), @(2), @(3)],
	[Z,  @(4), @(5).cond(e -> (e<=@(2).val) or (e >= @(3).val)) ]
      ],
 e -> [	II(@(1).val, (@(2).val-@(5).val) mod @(1).val,
           # this makes 0->size, because size mod size = 0
	   let(last:=(@(3).val-@(5).val) mod @(1).val,
	       When(last=0, @(1).val, last))) ]
);

# full_interval_diag * mat = mat
ARule( Compose,
      [ [Diag, [II, @(1), zero, @(2).cond(e->e=@(1).val)]], @(3) ],
 e -> [ @(3).val ]
);

# mat * full_interval_diag  = mat
ARule( Compose,
      [ @(3), [Diag, [II, @(1), zero, @(2).cond(e->e=@(1).val)]] ],
 e -> [ @(3).val ]
);

# ================================================================
# Sum of Monomials
# ================================================================
# *******************************************************************
Class( RulesGMonTensorPullOut, RuleSet );
# *******************************************************************

# GMon splitting
# GMon(fTensor, fTensor, fTensor) -> Tensor(GMon, GMon, ...)
#
tensor_ch := o -> When(IS_A(o, fTensor), o.children(), [o]);

Rule([GMon,
        @(1),
        @(2),
        @(3).cond(e ->
	    let(p := @(1).val, v := @(2).val, s := e,
	        ForAny([p,v,s], IS(fTensor)) and
	        compat_tensor_mons(tensor_ch(p), tensor_ch(v), tensor_ch(s))))],
 e -> let(
     cp := tensor_ch(@(1).val),
     cv := tensor_ch(@(2).val),
     cs := tensor_ch(@(3).val),
     Tensor( merge_tensor_mons(cp, cv, cs, GMon) ))
);

# ================================================================
# GMon -> SGMon
#  (delta * gath)
# ================================================================
# *******************************************************************
Class( RulesGMonToSGMon, RuleSet );
# *******************************************************************

Rule( [ GMon, @(1), @(2), full_interval(@(3),@(4)) ],
    e -> SGMon(fId(domain(@(1).val)), # scatter
	       @(1).val, # gather
	       @(2).val, # vdiag,
	       @(3).val));

## Diag(II) * Diag(..) * Gath(fId)
Rule( [ GMon, @(1,fId), @(2), @(3,II) ],
 e -> let(left := @(3).val.params[2],
          right:= @(3).val.params[3],
          size := @(3).val.params[1],
	  fadd := fAdd(size, right-left, left),
    SGMon(fadd, fadd,
          fCompose(@(2).val, fadd), II(right-left)))
);

## Diag(II) * Diag(..) * Gath(Z)
Rule( [ GMon,
	[Z, @n, @z],
	@(1),
	[II, @, @l, @(3).cond(r -> (@z.val <= @n.val-r) or
		                   (@z.val >= @n.val-@l.val)) ] ],
e -> let(size  := @n.val,
	 left  := @l.val,
         right := @(3).val,
	 shift := @z.val,
	 scat := fAdd(size, right-left, left),
	 gath := fAdd(size, right-left, (left + shift) mod size),
    SGMon(scat, gath,
	  fCompose(@(1).val, scat), II(domain(gath))))
);

## Diag(II o fBase) * Gath(fBase)
Rule([GMon,
	@(1,fBase),
	@(2),
	[fCompose, @(3,II), same_params(@(4),fBase,@(1)) ] ],

 e -> let(left := @(3).val.params[2], right := @(3).val.params[3],
          j    := @(1).val.params[2],
	  gath := @(1).val,
	  COND(leq(left, j, right-1),
	       SGMon(fId(1), gath, @(2).val, II(1)),
	       O(1, range(gath))))
);

## Diag(II o fBase) * Gath(fId(1))
Rule([skip,
	@(1,fId,e->range(e)=1),
	@(2),
	[fCompose, @(3,II), @(4,fBase) ] ],

 e -> let(left := @(3).val.params[2], right := @(3).val.params[3],
          j    := @(4).val.params[2],
	  COND(leq(left, j, right-1),
	       SGMon(fId(1), fId(1), @(2).val, II(1)),
	       O(1, 1)))
);

# NOTE:Make direct sum rules more general (... + p + ...) o fBase
#
## Diag(leftII o fBase) * Gath((p + ...) o fBase)
Rule([GMon,
        [fCompose,
	    [fDirsum, @(1).cond(e -> is_perm(e)), ...],
	    @(2,fBase) ],
	@(3),
	[fCompose, left_interval(@(10), @(20), @(4).cond(e->range(@(1).val)=e)),
	           same_params(@(5), fBase, @(2))] ],

 e -> let(right := @(4).val,
          size  := @(2).val.params[1],
          j     := @(2).val.params[2],
	  COND(leq(0, j, right-1),
	       SGMon(
		   fId(1),
		   fCompose(fAdd(size, right, 0), @(1).val, fBase(right, j)),
		   @(3).val,
		   II(1)),
	       O(1, size)))
);

## Diag(rightII o fBase) * Gath((... + p) o fBase)
Rule([GMon,
        [fCompose,
	    [fDirsum, ..., @(1).cond(e -> is_perm(e))],
	    @(2,fBase) ],
	@(3),
	[fCompose, right_interval(@(0), @(4), @(5).cond(e->range(@(1).val) = @(4).val - e)),
	           same_params(@(6), fBase, @(2)) ] ],

 e -> let(left := @(4).val,
          size  := @(2).val.params[1],
          j     := @(2).val.params[2],
	  COND(leq(left, j, size-1),
	       SGMon(
		   fId(1),
		   fCompose(fAdd(size, size-left, left), @(1).val, fBase(size-left, j-left)),
		   @(3).val,
		   II(1)),
	       O(1, size)))
);

# ================================================================
# Going back from Tensor(SGMon, ...) into SGMon(fTensor,...)
# ================================================================

Rule([@(0,Tensor), ..., @(1,COND), ...],
    e -> let( ch := e.children(),
	      left := Left(...),
	      right := Right(...),
	      leftch := ch{[1..left]},
	      rightch := ch{[right..Length(ch)]},
	      OP := ObjId(@(0).val),
	      COND(@(1).val.cond,
		  List(@(1).val.children(),
		       c -> OP(Concatenation(Copy(leftch), [c], Copy(rightch)))))));

# *******************************************************************
Class( RulesSGMonTensorPullIn, RuleSet );
# *******************************************************************

Rule( @(1,Tensor,CHILDREN_ARE(SGMon)),
 e -> let( ch:=e.children(),
     SGMon(
	 fTensor(List(ch, c->c.scat)),
	 fTensor(List(ch, c->c.gath)),
	 fTensor(List(ch, c->c.vdiag)),
	 fTensor(List(ch, c->c.sdiag)))));

# ================================================================
# Finalization
# ================================================================
Class( RulesFinalize, RuleSet );

#z:=SUMAcc(Gath(fId(4)), Scat(fAdd(4,3,0)) * Gath(fAdd(4,3,0)));

Elim := (elim_func, s) -> When(elim_func(s), [], s);

ARule(SUMAcc,
    [ [Gath, @R],
      [Compose,
	  [@(1,Scat), [@(2,fAdd), @, @, @k]],
	  [@(3,Gath), @S]] ],

    e -> let(n   := domain(@R.val), N   := range(@R.val),
	     nmk := domain(@S.val), nm  := nmk + @k.val,
	     m := n - nm,
    [ SUM( Elim(e->@k.val=0, Scat(fAdd(n, @k.val, 0)) * Gath(fCompose(@R.val, fAdd(n, @k.val, 0)))),
	   Scat(Copy(@(2).val)) *
	   SUMAcc(
	       Gath(fCompose(@R.val, fAdd(n, nmk, @k.val))),
	       Gath(@S.val)),
	   Elim(e->m=0, Scat(fAdd(n, m, nm)) * Gath(fCompose(@R.val, fAdd(n, m, nm))))
      ) ])
);

ARule(Compose,
    [ [Scat, @W],
      @(1, SUMAcc, e->ForAll(e.children(), c->ObjId(c)=Gath and domain(c.func)>1)) ],

  e -> [ let(n := domain(@W.val),
	     i := Ind(n),
	     fbase := fBase(n, i),
	  ISum(i, n,
	      Scat(fCompose(@W.val, fbase)) *
	      SUMAcc(List(@(1).val.children(), c -> Gath(fCompose(c.func, fbase)))))) ]
);

ARule(Compose,
    [ [Scat, @W], [Diag, @D], [Gath, @R] ],

  e -> [ let(n := domain(@R.val),
	     i := Ind(n),
	     fbase := fBase(n, i),
	  ISum(i, n,
	      Scat(fCompose(@W.val, fbase)) *
	      Blk1(@D.val.lambda().at(i)) *
	      Gath(fCompose(@R.val, fbase)))) ]
);

ARule(Compose, [@(1, O), @(2)], e->[ @(1).val ]);

Rule([Compose, [Scat, @(1)], @(2, O)], e->
    let(n := @(1).val.domain(),
	i := Ind(n),
	fbase := fBase(n,i),
	ISum(i, n,
	     Scat(fCompose(@(1).val, fbase)) *
	     Blk1(0) *
	     Gath(fCompose(fId(Cols(@(2).val)), fbase))))
);


text:
package main

import (
	"github.com/containous/yaegi/interp"
)

func main() {
	i := interp.New(interp.Opt{})
	i.Use(interp.ExportValue, interp.ExportType)
	i.Eval(`import "github.com/containous/yaegi/interp"`)
	i.Eval(`i := interp.New(interp.Opt{})`)
	i.Eval(`i.Eval("println(42)")`)
}

// Output:
// 42

ext: gi

max_stars_repo_name: blasrodri/yaegi

max_issues_repo_name: blasrodri/yaegi

max_forks_repo_name: blasrodri/yaegi


text:
# Built-in
UnorderedTuples(["iced", "jam", "plain"], 2);

ext: gap

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


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

VarFrontier := (e, depth) -> Cond(
    IsBound(e.isExpComposite) and e.isExpComposite,
        Union(List(VarArgsExp(e), x->VarFrontier(x, depth))),
    IsVar(e),
        When(depth=0,
            Set([e]),
            Union([e], Union(List(PredLoc(e), x->VarFrontier(x, depth-1))))),
    Set([])
);

IsXRef := exp -> Intersection([var.table.X,var.table.X1,var.table.XY1], VarFrontier(exp, -1)) <> [];

IsYRef := exp -> Intersection([var.table.Y,var.table.Y1,var.table.XY1], VarFrontier(exp, -1)) <> [];

CollectLoads := cmd -> Collect(cmd, [assign, @, @(1, [deref, nth], IsXRef)]);

CollectStores := cmd -> Collect(cmd, [assign, @(1, [deref, nth], IsYRef), @]);


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

# Option Records
# ==============
# MP, BS, from 01/25/01, GAPv3.4.4

# general default search options
SEARCH_DEFAULTS :=
    rec(
		timeLimit          := false,
		globalUnrolling    := false,
		globalUnrollingMin := 8,
		globalUnrollingMax := 64,
		timeBaseCases      := true
    );

# default DP search options
DP_DEFAULTS :=
    rec(
        nBest     := 1,
        optimize  := "minimize",
        verbosity := 3,
	    DPVec     := false,
	    DPVecVlen := 4
    );

#F Search Options Record
#F ---------------------

#F Search options records keep search options that are not specific to any
#F particular search algorithm.  These are options that cause the search
#F algorithms to search over different implementations instead of just how
#F the search algorithm progresses.
#F
#F These options get incorporated into a search algorithm's options record.
#F

#F PrintSpecSearchOptionsRecord()
#F   prints the specification for the general search options
#F   only to be called by specific search algorithm's PrintSpec
#F

PrintSpecSearchOptionsRecord := function()
   Print("  timeLimit := false | <minutes>,\n");
   Print("  globalUnrolling := true | false,\n");
   Print("  globalUnrollingMin := <positive int>,\n");
   Print("  globalUnrollingMax := <positive int>,\n");
end;

#F CheckSearchOptionsRecord( <search-options-record> )
#F   checks whether <search-options-record> is a valid search options record
#F

CheckSearchOptionsRecord := function( R )
   local r;

   if not IsRec(R) then
    Error("<R> must be an search options record");
  fi;

  # check fields
  for r in RecFields(R) do
     if r = "timeLimit" then
        if not ( R.(r) = false or ( IsInt(R.(r)) and R.(r) > 0 ) ) then
       Error( "timeLimit must be either false or a positive integer" );
    fi;

     elif r = "globalUnrolling" then
        if not IsBool(R.(r)) then
       Error( "Search option globalUnrolling must be true or false" );
    fi;
     elif r = "globalUnrollingMin" then
        if not IsInt(R.(r)) then
       Error( "globalUnrollingMin is not an integer" );
    fi;
     elif r = "globalUnrollingMax" then
        if not IsInt(R.(r)) then
       Error( "globalUnrollingMax is not an integer" );
    fi;
     elif r = "timeBaseCases" then
        if not IsBool(R.(r)) then
       Error( "Search option timeBaseCases must be true or false" );
    fi;
     fi;
  od;

  return true;
end;

#F MergeSearchOptionsRecord( <search-options-record> )
#F   merges <search-options-record> with the default values
#F

MergeSearchOptionsRecord := function( R )
   local OR, r;

   CheckSearchOptionsRecord(R);

   OR := ShallowCopy(R);
   for r in RecFields(SEARCH_DEFAULTS) do
      if not IsBound( OR.(r) ) then
         OR.(r) := SEARCH_DEFAULTS.(r);
      fi;
   od;

   if OR.globalUnrolling = true and
      OR.globalUnrollingMin > OR.globalUnrollingMax then
      Error( "globalUnrollingMin > globalUnrollingMax" );
   fi;

   return OR;
end;

#F DP Options Record
#F -----------------

#F The DP Options Records maintain options for the Dynamic Programming search
#F algorithm.  This is a combination of both search options specific to DP
#F as well as general search options.

#F PrintSpecDPOptionsRecord()
#F   prints the specification for the DP search options
#F

PrintSpecDPOptionsRecord := function()
   Print("rec(\n");
   PrintSpecSearchOptionsRecord();
   Print("  nBest := <positive integer>,\n");
   Print("  optimize := \"minimize\" | \"maximize\",\n");
   Print("  hashTable := <hashTable>,\n");
   Print("  verbosity := <non-negative integer>\n");
   Print(");\n");
end;

#F CheckDPOptionsRecord( <DP-options-record> )
#F   checks to see if DP-options-record is valid
#F

CheckDPOptionsRecord := function( R )
   local r;

   CheckSearchOptionsRecord(R);
   for r in RecFields(R) do
      if r = "nBest" then
         if not ( IsInt(R.(r)) and R.(r) > 0 ) then
        Error( "nBest must be a positive integer" );
     fi;
      elif r = "verbosity" then
         if not ( IsInt(R.(r)) and R.(r) >= 0 ) then
        Error( "verbosity must be a non-negative integer" );
     fi;
      elif r = "optimize" then
         if not R.(r) in ["minimize","maximize"] then
        Error( "optimize must be to minimize or maximize" );
     fi;
      elif r = "hashTable" then
         if not IsHashTable( R.(r) ) then
        Error( "hashTable must be a valid hashTable" );
         fi;
      elif r = "hashTableBases" then
         if not IsHashTable( R.(r) ) then
        Error( "hashTableBases must be a valid hashTable" );
         fi;
      elif r = "breakdownRules" then
         if not IsRec( R.(r) ) then
        Error( "breakdownRules must be a record" );
         fi;
      elif not r in RecFields(SEARCH_DEFAULTS) and  r <> "DPVec" and r <> "measureFunction" and r <> "wrap"
      and not IsSystemRecField(r) then
         Error( "Unknown DPOptionsRecord field <r>" );
      fi;
   od;
end;

#F MergeDPOptionsRecord( <DP-options-record> )
#F   merges <DP-options-record> with the defaults
#F

MergeDPOptionsRecord := function( R )
   local OR, r;

   CheckDPOptionsRecord(R);
   OR := MergeSearchOptionsRecord(R);
   for r in RecFields(DP_DEFAULTS) do
      if not IsBound( OR.(r) ) then
         OR.(r) := DP_DEFAULTS.(r);
      fi;
   od;

   return OR;
end;

# Implement Options Record
# ------------------------

#F The Implement Options Records maintain options for timed search.
#F

# Utility for Search algorithms to check dataType in SPLOpts
# ----------------------------------------------------------

#F SearchCheckDataType( <spl>, <SPL-options-record> )
#F   Checks to see that the dataType is set properly for the given spl.
#F   Sets it if no default was given.
#F

SearchCheckDataType := function( spl, SPLOpts )
   if not IsBound( SPLOpts.dataType ) then
      Error( "SPLOpts.dataType not bound" );
   fi;

   if SPLOpts.dataType = "no default" then
      if IsRealSPL(spl) then
         SPLOpts.dataType := "real";
      else
         SPLOpts.dataType := "complex";
      fi;
   elif SPLOpts.dataType in ["real","complex"] then
       ;
      # YEVGEN: This should not be too ambitious with error checking
      # we currently use real datatype for complex transforms (RC operator is applied
      # to get a real formula).

      #if (not IsRealSPL(spl)) and SPLOpts.dataType = "real" then
      #   Error("Real data type specified but spl is complex");
      #fi;
   else
      Error("SPLOpts.dataType is not \"real\" or \"complex\"");
   fi;
end;


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

NewRulesFor(TTensorInd, rec(
#   base cases
#   I x A
    dsA_base_vec_push := rec(
        info := "IxA base",
        forTransposition := false,
        applicable := nt -> IsParPar(nt.params) and nt.isTag(1, AVecReg),
        children := nt -> let(_krnl := nt.params[1].withTags(nt.getTags()), krnl := When(_krnl.isReal(), _krnl.setWrap(VWrapId), _krnl.setWrap(VWrapTRC(nt.firstTag().isa))),
            [[ krnl, InfoNt(nt.params[2]) ]]),
        apply := (nt, c, cnt) -> IDirSum(cnt[2].params[1], c[1])
    ),
#   A x Iv
    L_dsA_L_base_vec := rec(
        info := "AxI base",
        forTransposition := false,
        applicable := nt -> IsVecVec(nt.params) and nt.isTag(1, AVecReg) and nt.params[2].range = nt.getTags()[1].v,
        children := nt -> let(jv := Ind(nt.getTags()[1].v),
            [[ SubstVars(Copy(nt.params[1]), rec((nt.params[2].id) := jv)).setWrap(VWrap(nt.firstTag().isa)).withTags(Drop(nt.getTags(), 1)),
            InfoNt(jv) ]]),
        apply := (nt, c, cnt) -> let(v := nt.getTags()[1].v, jv := cnt[2].params[1],
            VTensorInd(c[1], jv))
    ),
#   A x In
    L_dsA_L_vec := rec(
        info := "AxI base",
        forTransposition := false,
        applicable := nt -> IsVecVec(nt.params) and nt.isTag(1, AVecReg) and nt.params[2].range > nt.getTags()[1].v,
        children := nt -> let(v := nt.getTags()[1].v, jv := Ind(nt.getTags()[1].v), j := Ind(nt.params[2].range/nt.getTags()[1].v),
            [[ TTensorInd(
                SubstVars(Copy(nt.params[1]), rec((nt.params[2].id) := j*V(v)+jv)),
                j, AVec, AVec).setWrap(VWrap(nt.firstTag().isa)).withTags(Drop(nt.getTags(), 1)), InfoNt(jv) ]]),
        apply := (nt, c, cnt) -> let(jv := cnt[2].params[1],
            VTensorInd(c[1], jv))
    )
));


text:
# Spiral matrix with numbers 1 .. n2, more natural in GAP
SpiralMatrix := function(n)
  local i, j, k, di, dj, p, vi, vj, imin, imax, jmin, jmax;
  a := NullMat(n, n);
  vi := [ 1, 0, -1, 0 ];
  vj := [ 0, 1, 0, -1 ];
  imin := 0;
  imax := n;
  jmin := 1;
  jmax := n + 1;
  p := 1;
  di := vi[p];
  dj := vj[p];
  i := 1;
  j := 1;
  for k in [1 .. n*n] do
    a[j][i] := k;
    i := i + di;
    j := j + dj;
    if i < imin or i > imax or j < jmin or j > jmax then
      i := i - di;
      j := j - dj;
      p := RemInt(p, 4) + 1;
      di := vi[p];
      dj := vj[p];
      i := i + di;
      j := j + dj;
      if p = 1 then
        imax := imax - 1;
      elif p = 2 then
        jmax := jmax - 1;
      elif p = 3 then
        imin := imin + 1;
      else
        jmin := jmin + 1;
      fi;
    fi;
  od;
  return a;
end;

PrintArray(SpiralMatrix(5));
# [ [   1,   2,   3,   4,   5 ],
#   [  16,  17,  18,  19,   6 ],
#   [  15,  24,  25,  20,   7 ],
#   [  14,  23,  22,  21,   8 ],
#   [  13,  12,  11,  10,   9 ] ]


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

Declare(VContainer);

###############################################################################
#F VContainer(<spl>, <isa>)
#F
#F   self.child(1) -- spl
#F   self.isa      -- isa
#F
Class(VContainer, VRC, rec(
    #-----------------------------------------------------------------------
    new := (self, spl, isa) >> SPL(WithBases(self, rec(
        _children := [spl], isa := isa))).setDims(),
    #-----------------------------------------------------------------------
    dims := self >> self.child(1).dims(),
    #-----------------------------------------------------------------------
    toAMat := (self) >> AMatMat(MatSPL(self.child(1))),
    #-----------------------------------------------------------------------
    isReal := self >> self.child(1).isReal(),
    #-----------------------------------------------------------------------
    conjTranspose := self >> ObjId(self)(self.child(1).conjTranspose(), self.isa),
    #-----------------------------------------------------------------------
    # NOTE: make sure rChildren properly exposes all parameters
    from_rChildren := (self, rch) >> VContainer(rch[1], self.isa),
    #-----------------------------------------------------------------------
    print := (self, i, is) >> Print(self.name, "(\n", Blanks(i+is),
	self.child(1).print(i+is,is), ",\n", Blanks(i+is), self.isa, ")",
	self.printA()),
    #-----------------------------------------------------------------------
    unroll := self >> self,
    #-----------------------------------------------------------------------
    transpose := self >> VContainer(self.child(1).transpose(), self.isa),
    #-----------------------------------------------------------------------
    vcost := self >>  self.child(1).vcost()
));

Class(VirtualPad, BaseContainer, SumsBase, rec(
    #-----------------------------------------------------------------------
    new := (self, rows, cols, ch) >> SPL(WithBases(self, rec(
        dimensions := [rows, cols], _children := [ch]))).setDims(),
    #-----------------------------------------------------------------------
    dims := self >> self.dimensions,
    #-----------------------------------------------------------------------
    toAMat := self >> let(m := MatSPL(self.child(1)), ch := self.child(1),
	AMatMat(
	    List(m, x -> x :: Replicate(Cols(self)-Cols(ch), 0)) ::
	    Replicate(Rows(self)-Rows(ch), Replicate(Cols(self), 0)))),
    #-----------------------------------------------------------------------
    rChildren := self >> [self.dimensions[1], self.dimensions[2], self._children[1]],
    #-----------------------------------------------------------------------
    rSetChild := meth(self, n, what)
        if   n=1 then self.dimensions[1] := Checked(IsIntSym(what), what);
        elif n=2 then self.dimensions[2] := Checked(IsIntSym(what), what);
	elif n=3 then self._children[1]  := Checked(IsSPL(what), what);
	else Error("<n> must be in [1..3]");
	fi;
    end,
));


text:
# Dirty implementation
# Only P3 format, an image is a list of 3 matrices (r, g, b)
# Max color is always 255
WriteImage := function(name, img)
  local f, r, g, b, i, j, maxcolor, nrow, ncol, dim;
  f := OutputTextFile(name, false);
  r := img[1];
  g := img[2];
  b := img[3];
  dim := DimensionsMat(r);
  nrow := dim[1];
  ncol := dim[2];
  maxcolor := 255;
  WriteLine(f, "P3");
  WriteLine(f, Concatenation(String(ncol), " ", String(nrow), " ", String(maxcolor)));
  for i in [1 .. nrow] do
    for j in [1 .. ncol] do
      WriteLine(f, Concatenation(String(r[i][j]), " ", String(g[i][j]), " ", String(b[i][j])));
    od;
  od;
  CloseStream(f);
end;

PutPixel := function(img, i, j, color)
  img[1][i][j] := color[1];
  img[2][i][j] := color[2];
  img[3][i][j] := color[3];
end;

GetPixel := function(img, i, j)
  return [img[1][i][j], img[2][i][j], img[3][i][j]];
end;

NewImage := function(nrow, ncol, color)
  local r, g, b;
  r := color[1] + NullMat(nrow, ncol);
  g := color[2] + NullMat(nrow, ncol);
  b := color[3] + NullMat(nrow, ncol);
  return [r, g, b];
end;

# Reproducing the example from Wikipedia
black := [ 0, 0, 0 ];
g := NewImage(2, 3, black);
PutPixel(g, 1, 1, [255, 0, 0]);
PutPixel(g, 1, 2, [0, 255, 0]);
PutPixel(g, 1, 3, [0, 0, 255]);
PutPixel(g, 2, 1, [255, 255, 0]);
PutPixel(g, 2, 2, [255, 255, 255]);
PutPixel(g, 2, 3, [0, 0, 0]);
WriteImage("example.ppm", g);


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

###############################################################################
##
#R  IsSelfSimFamilyRep
##
##  Any object from category SelfSimFamily which is created here, lies in a
##  category IsSelfSimFamilyRep.
##  Family of SelfSim object  < a >  contains all essential information about
##  (mathematical) automaton which generates group containing  < a > :
##  it contains automaton, properties of automaton and group generated by it,
##  etc.
##  Also family contains group generated by states of underlying automaton.
##
DeclareRepresentation("IsSelfSimFamilyRep",
                      IsComponentObjectRep and IsAttributeStoringRep,
                      [ "freegroup",      # IsSelfSim objects store words from this group
                        "freegens",       # list [f1, f2, ..., fn, f1^-1, ..., fn^-1, 1]
                                          # where f1..fn are generators of freegroup,
                                          # n is numstates; 1 is stored if and only if
                                          # trivstate is not zero.
                                          # Some fi^-1 may be missing if corresponding
                                          # generator is not invertible (but the list still
                                          # has the length of 2n + 1).
                        "recurgens",      # Generators of the group, list of length n.
                        "numstates",      # number of non - trivial generating states
                        "deg",
                        "trivstate",      # 0 or 2*numstates + 1
                        "isgroup",        # whether all generators are invertible
                        "names",          # list of non - trivial generating states
                        "recurlist",      # the automaton table, states correspond to freegens
                        "oldstates",      # mapping from states in the original table used to
                                          # define the group to the states in recurlist:
                                          # SelfSim(fam!.freegens[oldtstates[k]], fam) is the element
                                          # which corresponds to k-th state in the original automaton
                        "rws",            # rewriting system
                        "use_rws",        # whether to use rewriting system in multiplication
                        "use_contraction" # whether to use contraction in IsOne
                      ]);

BindGlobal("AG_FixRecurList",
function(list, oldstates, names)
  local i, j, k, reduced, deg, isgroup, numstates, perm, trivstate, trivstates, new_trivstates;

  deg := Length(list[1]) - 1;
  isgroup := true;
  trivstate := 0;
  trivstates := [];
  numstates := Length(list);

  # convert list to a "canonical" form

  for i in [1..numstates] do
    for j in [1..deg] do
      if not IsList(list[i][j]) then list[i][j] := [list[i][j]]; fi;
    od;
  od;

  # find if there are any "obviously" trivial states

  repeat
    new_trivstates := [];
    #  find new trivial states (new ones may appear in every iteration)
    for i in [1..Length(list)] do
      if (not i in trivstates) and AG_IsObviouslyTrivialStateInList(i, list) then
        Add(new_trivstates, i);
        numstates := numstates - 1;
      fi;
    od;

    for trivstate in new_trivstates do
      for i in [1..Length(list)] do
        for j in [1..deg] do
      #  remove trivstate from everywhere
          list[i][j] := Filtered(list[i][j], x -> AbsInt(x) <> trivstate);
      #  freely reduce words
          repeat
            reduced := true;
            for k in [1..Length(list[i][j]) - 1] do
              if list[i][j][k] = -list[i][j][k + 1] then
                Remove(list[i][j], k);
                Remove(list[i][j], k);
                reduced := false;
                break;
              fi;
            od;
          until reduced;
        od;
      od;
    od;
    Append(trivstates, new_trivstates);
  until new_trivstates = [];

  # move trivial state to the end
  Sort(trivstates);
  trivstates := Reversed(trivstates);

  for trivstate in trivstates do
    for i in [1..Length(list)] do
      for j in [1..deg] do
    #  remove trivstate from everywhere
        Apply(list[i][j], function(x)
                            if x > trivstate then return x - 1;
                            elif x < -trivstate then return x + 1;
                            fi;
                            return x;
                          end);
      od;
    od;

    for i in [1..Length(oldstates)] do
      if oldstates[i] = trivstate then
        oldstates[i] := 2*numstates + 1;
      elif oldstates[i] > trivstate then
        oldstates[i] := oldstates[i] - 1;
      fi;
    od;
    Remove(list, trivstate);
    Remove(names, trivstate);
  od;

  if trivstates <> [] then
    trivstate := 2*numstates + 1;
    list[trivstate] := List([1..deg], k -> []);
    list[trivstate][deg + 1] := ();
  fi;

  # add inverses of states
  for i in [1..numstates] do
    if AG_IsInvertibleStateInList(i, list) then
      list[i + numstates] := [];

      list[i][deg + 1] := AG_PermFromTransformation(list[i][deg + 1]);

      perm := list[i][deg + 1];
      list[i + numstates][deg + 1] := perm^-1;
      for j in [1..deg] do
        list[i + numstates][j] := -Reversed(list[i][j^(perm^-1)]);
      od;
    else
      isgroup := false;
      if list[i][deg + 1]^-1<>fail then
        list[i][deg + 1] := AG_PermFromTransformation(list[i][deg + 1]);
      fi;
    fi;
  od;

  return [list, numstates, trivstate, isgroup];
end);

###############################################################################
##
#M  SelfSimFamily( <list>, <names>, <bind_vars> )
##
InstallMethod(SelfSimFamily, "for [IsList, IsList, IsBool]",
              [IsList, IsList, IsBool],
function (list, names, bind_global)
  local deg, tmp, trivstate, numstates, i,
        freegroup, freegens, a, family, oldstates,
        isgroup;

  if not AG_IsCorrectRecurList(list, false) then
    Error("in SelfSimFamily(IsList, IsList, IsString):\n",
          "  given list is not a correct list representing self-similar group\n");
  fi;

# 1. make a local copy of arguments, since they will be modified and put into the result

  list := StructuralCopy(list);
  names := StructuralCopy(names);
  deg := Length(list[1]) - 1;

# 2. Find trivial states, permute states

  oldstates := [1..Length(list)];
  tmp := AG_FixRecurList(list, oldstates, names);
  list := tmp[1];
  numstates := tmp[2];
  trivstate := tmp[3];
  isgroup := tmp[4];

# 3. Create FreeGroup and FreeGens

  freegroup := FreeGroup(names);
  freegens := ShallowCopy(FreeGeneratorsOfFpGroup(freegroup));
  for i in [1..numstates] do
    if IsBound(list[i + numstates]) then
      freegens[i + numstates] := freegens[i]^-1;
    fi;
  od;
  if trivstate <> 0 then
    freegens[trivstate] := One(freegroup);
  fi;

# 4. Create family

  if isgroup then
    family := NewFamily("SelfSimFamily",
                        IsInvertibleSelfSim,
                        IsInvertibleSelfSim,
                        IsSelfSimFamily and IsSelfSimFamilyRep);
  else
    family := NewFamily("SelfSimFamily",
                        IsSelfSim,
                        IsSelfSim,
                        IsSelfSimFamily and IsSelfSimFamilyRep);
  fi;

  family!.isgroup := isgroup;
  family!.deg := deg;
  family!.numstates := numstates;
  family!.trivstate := trivstate;
  family!.names := names;
  family!.freegroup := freegroup;
  family!.freegens := freegens;
  family!.recurlist := list;
  family!.oldstates := oldstates;
  family!.use_rws := false;
  family!.rws := fail;
  family!.use_contraction := false;

  SetIsActingOnBinaryTree(family, deg = 2);
  SetDegreeOfTree(family, deg);
  SetTopDegreeOfTree(family, deg);

  family!.recurgens := [];
  for i in [1..Length(list)] do
    if IsBound(list[i]) then
      family!.recurgens[i]  :=
        __AG_CreateSelfSim(family, freegens[i],
                        List([1..deg], j -> AssocWordByLetterRep(FamilyObj(freegens[1]), list[i][j])),
                        list[i][deg + 1],
                        i > numstates or IsBound(list[i + numstates]));
    fi;
  od;

  # XXX It's evil to bind global names, consider AssignGeneratorVariables
  # XXX Check whether names are actually valid names for variables
  if bind_global then
    for i in [1..family!.numstates] do
      if IsBoundGlobal(family!.names[i]) then
        UnbindGlobal(family!.names[i]);
      fi;
      BindGlobal(family!.names[i], family!.recurgens[i]);
      MakeReadWriteGlobal(family!.names[i]);
    od;
    #BindGlobal(AG_Globals.identity_symbol, One(family));
    #MakeReadWriteGlobal(AG_Globals.identity_symbol);
  fi;

  return family;
end);

###############################################################################
##
#M  SelfSimFamily( <list> )
##
InstallMethod(SelfSimFamily, "for [IsList]", [IsList],
function(list)
  return SelfSimFamily(list, false);
end);

###############################################################################
##
#M  SelfSimFamily( <list>, <bind_vars> )
##
InstallMethod(SelfSimFamily, "for [IsList, IsBool]", [IsList, IsBool],
function(list, bind_vars)
  if not AG_IsCorrectRecurList(list, false) then
    Error("in SelfSimFamily(IsList):\n",
          "  given list is not a correct list representing self-similar group\n");
  fi;

  return SelfSimFamily(list,
                     List([1..Length(list)],
                          i -> Concatenation(AG_Globals.state_symbol, String(i))),
                     bind_vars);
end);

###############################################################################
##
#M  SelfSimFamily( <list>, <names> )
##
InstallMethod(SelfSimFamily, "for [IsList, IsList]", [IsList, IsList],
function(list, names)
  if not AG_IsCorrectRecurList(list, false) then
    Error("in SelfSimFamily(IsList, IsList):\n",
          "  given list is not a correct list representing automaton\n");
  fi;

  return SelfSimFamily(list, names, AG_Globals.bind_vars_autom_family);
end);

###############################################################################
##
#M  One( <fam> )
##
InstallMethod(One, "for [IsSelfSimFamily]", [IsSelfSimFamily],
function(fam)
  return __AG_CreateSelfSim(fam, One(fam!.freegroup),
                         List([1..fam!.deg], i -> One(fam!.freegroup)),
                         (), true);
end);

###############################################################################
##
#M  GroupOfSelfSimFamily( <fam> )
##
InstallMethod(GroupOfSelfSimFamily, "for [IsSelfSimFamily]",
              [IsSelfSimFamily],
function(fam)
  local g;

  if not fam!.isgroup then
    return fail;
  fi;

  if fam!.numstates > 0 then
    g := GroupWithGenerators(fam!.recurgens{[1..fam!.numstates]});
  else
    g := Group(One(fam));
  fi;

  SetUnderlyingSelfSimFamily(g, fam);
  SetIsGroupOfSelfSimFamily(g, true);
  # XXX
  SetGeneratingRecurList(g, fam!.recurlist);
  SetRecurList(g, fam!.recurlist);
  SetDegreeOfTree(g, fam!.deg);
  SetTopDegreeOfTree(g, fam!.deg);
  SetIsActingOnBinaryTree(g, fam!.deg = 2);

  return g;
end);

###############################################################################
##
#M  SemigroupOfSelfSimFamily( <fam> )
##
InstallMethod(SemigroupOfSelfSimFamily, "for [IsSelfSimFamily]",
              [IsSelfSimFamily],
function(fam)
  local g;

  if fam!.trivstate <> 0 then
    if fam!.numstates = 0 then
      g := SemigroupByGenerators([One(fam!.recurgens[fam!.trivstate])]);
    else
      g := MonoidByGenerators(fam!.recurgens{[1..fam!.numstates]});
    fi;
  else
    g := SemigroupByGenerators(fam!.recurgens{[1..fam!.numstates]});
  fi;

  SetUnderlyingSelfSimFamily(g, fam);

  # XXX
  SetGeneratingRecurList(g, fam!.recurlist);
  SetRecurList(g, fam!.recurlist);

  SetDegreeOfTree(g, fam!.deg);
  SetTopDegreeOfTree(g, fam!.deg);
  SetIsActingOnBinaryTree(g, fam!.deg = 2);
  SetIsSelfSimilarSemigroup(g, true);

  return g;
end);

###############################################################################
##
#M  UnderlyingFreeMonoid( <G> )
##
InstallMethod(UnderlyingFreeMonoid, "for [IsSelfSimFamily]",
              [IsSelfSimFamily],
function(fam)
  local monoid;

  if fam!.numstates <> 0 then
    monoid := MonoidByGenerators(GeneratorsOfGroup(fam!.freegroup));
    SetSize(monoid, infinity);
  else
    monoid := MonoidByGenerators(fam!.freegens[1]);
    SetSize(monoid, 1);
  fi;

  return monoid;
end);

###############################################################################
##
#M  UnderlyingFreeGroup( <G> )
##
InstallMethod(UnderlyingFreeGroup, "for [IsSelfSimFamily]",
              [IsSelfSimFamily],
function(fam)
  return fam!.freegroup;
end);

###############################################################################
##
#M  IsObviouslyFiniteState( <G> )
##
##  returns `true' if there are no words longer than 1 in the wreath recursion
InstallMethod(IsObviouslyFiniteState, "for [IsSelfSimFamily]",
              [IsSelfSimFamily],
function(fam)
  local list, g, i;
  list := fam!.recurlist;
  for g in list do
    for i in [1..fam!.deg] do
      if Length(g[i]) > 1 then return false; fi;
    od;
  od;
  IsFiniteState(GroupOfSelfSimFamily(fam));
  return true;
end);

#############################################################################
##
#M  GeneratorsOfOrderTwo( <fam> )
##
InstallMethod(GeneratorsOfOrderTwo, "for [IsSelfSimFamily]", [IsSelfSimFamily],
function(fam)
  if not fam!.isgroup then
    Error("not all generators of the family are invertible");
  fi;
  return Filtered(GeneratorsOfGroup(GroupOfSelfSimFamily(fam)), g -> IsOne(g^2));
end);

###############################################################################
##
##  AG_AbelImagesGenerators(<fam>)
##
InstallMethod(AG_AbelImagesGenerators, "for [IsAutomFamily]",
              [IsSelfSimFamily],
function(fam)
  if not fam!.isgroup then
    Error("the underlying group is not generated by automorphisms");
  fi;
  return AG_AbelImageAutomatonInList(fam!.recurlist);
end);

#E

max_stars_repo_name: gap-packages/automgrp

max_issues_repo_name: gap-packages/automgrp

max_forks_repo_name: gap-packages/automgrp


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

neg.doPeel := true;

#F CopyPropagate(<code>)
#F     Performs copy propagation and strength reduction
#F     If <code> has IFs it MUST be in SSA form, i.e.,
#F     SSA(<code>) must be run first.
#F
Class(SubstVarsRules, RuleSet, rec(
    create := (self, varmap) >> Inherit(self,
	rec(
            rules := rec(
		v := Rule(@(1,var,e->IsBound(varmap.(e.id))), e -> varmap.(e.id))
	    )
	)).compile()
));

# apply linear normalization to index expressions, this will simplify cases like
# 3x + 2x = 5x
#
Class(SimpIndicesRules, RuleSet, rec(
    create := (self, opts) >> Inherit(self,
	rec(
	    rules := rec(
		simp_indices := Rule(@(1, opts.simpIndicesInside), x -> GroupSummandsExp(x))
	    )
	)).compile()
));

# Transforms nth(x, idx) to deref(x + idx)
Class(RulesDerefNth, RuleSet);
RewriteRules(RulesDerefNth, rec(
    deref_nth := Rule(
	[nth, @(1).cond( x -> not (x _is [Value, param]) or not IsBound(x.value)), @(2)],
	e -> let(
	    b := @(1).val, idx := @(2).val,
	    Cond(
		ObjId(idx) = add, deref(ApplyFunc(add, [b] :: idx.args)),
		ObjId(idx) = sub, deref(ApplyFunc(add, [b] :: [idx.args[1], neg(idx.args[2])])),
                                  deref(b + idx))))
));

# sorts the incoming array of assignments by the offset of the
# variable being assigned TO, rather than the one assigned FROM.
# improves locality hence cache performance.
_sortByIdx := function(array)
    local newarray;
    newarray := Copy(array);

    SortParallel(
        List(array, e -> Double(SubString(e.args[2].id, 2))),
        newarray
    );

    return newarray;
end;

CopyTab := function (t)
   local r, a;
   r := tab();
   for a in NSFields(t) do
      r.(a):=Copy(t.(a));
   od;
   return r;
end;

Class(CopyPropagate, CSE, rec(
    propagate := (self, cloc, cexp) >> let(cebase := ObjId(cexp),
        ((cebase in [var, Value, noneExp])
            or (self.propagateNth and cebase in [nth, deref] and IsValue(cexp.idx))
            or (IsBound(self.opts.autoinline) and IsBound(cloc.succ) and Length(cloc.succ)<=1)
            or (IsBound(cloc.succ) and Length(cloc.succ)=0))),

    # sreduce_and_subst(c.exp) can lead to inf loop - why?
    procRHS := (self, x) >> self.sreduce(self.subst3(self.sreduce(self.subst2(x)))),

    procLHS := (self, x) >> self.sreduce3(self.sreduce(self.subst2(x))),

    # NOTE: explain this
    prepVarMap := meth(self, initial_map, do_sreduce, do_init_scalarized, do_idx)
        local v, varmap, rs_subst, rs2, rs_deref, rs_idx;
        self.varmap := initial_map;

        if do_sreduce then
	    # do_idx is done initially, at the same time as scalarization
	    # deref prevents scalarization, so  must be disabled at that moment
	    rs_deref := When(not do_idx and self.opts.useDeref, RulesDerefNth, EmptyRuleSet);
	    rs_idx   := When(do_idx, SimpIndicesRules.create(self.opts), EmptyRuleSet);

            rs_subst := SubstVarsRules.create(self.varmap);
	    rs2 := MergedRuleSet(rs_subst, rs_deref);
	    rs_subst.__avoid__ := [Value];
	    rs_idx.__avoid__ := [Value];
	    rs2.__avoid__ := [Value];

            self.subst   := x -> SubstVars(x, self.varmap);
            self.subst2  := x -> BU(x, rs2);
            self.subst3  := x -> SubstBottomUpRules(BU(x, rs_subst), rs_idx.rules);
            self.sreduce := x -> SReduce(x, self.opts);
            self.sreduce3 := x -> SubstBottomUpRules(x, rs_idx.rules);
            self.sreduce_and_subst := MergedRuleSet(rs_subst, rs_deref, RulesStrengthReduce);
        else
            self.subst   := x -> SubstVars(x, self.varmap);
            self.subst2  := x -> SubstVars(x, self.varmap);
            self.subst3  := x -> SubstVars(x, self.varmap);
            self.sreduce := x -> x;
            self.sreduce_and_subst := self.subst;
        fi;

        if do_init_scalarized then
            for v in compiler.Compile.scalarized do
                self.varmap.(v.id) := noneExp(v.t);
            od;
        fi;
        return self;
    end,

    assumeSSA := false,
    afterSSA := self >> WithBases(self, rec(assumeSSA:=true)),

    closeVarMap := meth(self, other, newcmds)
        local v;
        for v in UserNSFields(other.varmap) do
            if (not IsBound(self.varmap.(v)) or self.varmap.(v) <> other.varmap.(v)) and SuccLoc(var(v))<>[] then
                if self.assumeSSA then ;
                #    self.varmap.(v) := other.subst(other.varmap.(v));
                #    PrintLine("close ", v, " => ", other.subst(other.varmap.(v)));
                else
                    Unbind(self.varmap.(v));
                    Add(newcmds.cmds, assign(var(v), other.subst(other.varmap.(v))));
                fi;
            fi;
        od;
        #Print("-----\n");
    end,

    init := meth(self, opts)
        self.opts := opts;
        self.prepVarMap(tab(), true, true, false);
        self.doScalarReplacement := opts.doScalarReplacement;
        self.propagateNth := opts.propagateNth;
        self.flush();
        return self;
    end,

    initial := meth(self, code, opts)
        self.opts := opts;
        self.prepVarMap(tab(), true, true, true);
        self.doScalarReplacement := opts.doScalarReplacement;
        self.propagateNth := opts.propagateNth;
        self.flush();
        return self.copyProp(code);
    end,

    __call__ := (self, code, opts) >> self.init(opts).copyProp(code),

    flush := meth(self)
         self.csetab := tab();
         if (self.doScalarReplacement) then self.lhsCSE := CSE.init(); fi;
    end,

    fast := meth(self, code, opts)
        self.prepVarMap(tab(), false, false, false);
        self.doScalarReplacement := opts.doScalarReplacement;
        self.propagateNth := opts.propagateNth;
        self.opts := opts;
        self.flush();
        return self.copyProp(code);
    end,

    procAssign := meth(self, c, newcmds, live_out)
        local newloc, newexp, cid, cloc, op, varmap, lkup;
        varmap := self.varmap;

        # NOTE: generalize this, use .in()/.out()/.inout() somehow
        if IsBound(c.p) then
            cid := (loc, exp) -> ObjId(c)(loc, exp, c.p);
        else cid := ObjId(c); fi;

        # to my current understanding, if we assign a phi function, nothing can be propagated.
        # If you propagate inside a phi, you lose the branch information.
        if (ObjId(c.exp)=phi) then Add(newcmds, c); return; fi;

        # Run strength reduction, variable remapping, and all other rewrite rules
        newexp := self.procRHS(c.exp);
        cloc := When(not IsVar(c.loc) or (c.loc in c.op_inout()),
	             self.procLHS(c.loc), c.loc);

        # check if newexp was already computed
        lkup := self.cseLookup(newexp);
        if lkup <> false then newexp := lkup; fi;

	# if cloc is marked as live_out, save it in the list, so that its not kicked out
        When(IsVar(cloc) and IsBound(cloc.live_out), AddSet(live_out, cloc));

        # invalidate the LHS in the CSE tables
        if (self.cseLookup(cloc) <> false)
            then self.cseInvalidate(cloc); fi;
        if (self.doScalarReplacement and self.lhsCSE.cseLookup(cloc) <> false)
            then self.lhsCSE.cseInvalidate(cloc); fi;

        # propagate
        if IsVar(cloc) and cid=assign and self.propagate(cloc,newexp) then
            # this should not happen due to sreduce/subst above
            When(IsVar(newexp) and IsBound(varmap.(newexp.id)), Error("Should not happen"));
            varmap.(cloc.id) := newexp;

        # do not propagate
        else
            # NOTE: YSV: is this right???
            #        it seems that its correct, above case catches var=noneExp,
            #        so this looks like a mem-store, e.g.,  nth(..) = noneExp
            if ObjId(newexp)=noneExp then return; fi;

            # do not propagate AND cloc is a variable
            if IsVar(cloc) and not (cloc in c.op_inout()) then
                # propagate 'neg' like unary operators outwards
                if IsBound(newexp.doPeel) and newexp.doPeel then
                    op := ObjId(newexp);
                    newexp := newexp.args[1];
                    if self.propagate(cloc, newexp) then
                        if IsVar(newexp) and IsBound(varmap.(newexp.id)) then
                            varmap.(cloc.id) := op(varmap.(newexp.id));
                        else
                            varmap.(cloc.id) := op(newexp);
                        fi;
                    else
                        newloc := cloc.clone();
                        varmap.(cloc.id) := op(newloc);
                        Add(newcmds, cid(newloc, newexp));
                        # careful! cid can be any storeop, not always <assign>
                        When(cid=assign, self.cseAdd(newloc, newexp));
                    fi;
                else
                    # variable that is not propagated is remapped to a fresh name (to get code in SSA form)
                    newloc := cloc.clone();
                    varmap.(cloc.id) := newloc;
                    Add(newcmds, cid(newloc, newexp));   # cid == assign | assign_nop | ...
                    # careful! cid can be any storeop, not always <assign>
                    When(cid=assign, self.cseAdd(newloc, newexp));
                fi;

            #  do not propagate AND cloc is *not* a variable (ie. nth(X, i)) or it is an inout variable
            else
                if self.doScalarReplacement or not (self.propagateNth or IsVar(newexp) or IsValue(newexp)) then
                    # for non-variable newexp a fresh temporary var sXX is created to hold the result
                    newloc := var.fresh_t("s", newexp.t);
                    self.cseAdd(newloc, newexp);
                    Add(newcmds, assign(newloc, newexp));
                    newexp := newloc;
                fi;

                if self.doScalarReplacement then
                    When(cid<>assign or not (ObjId(cloc) in [nth,deref]),
                        Error("Scalar replacement can't handle <cloc> of type ", ObjId(cloc)));
                    self.cseAdd(newexp, cloc);
                    self.lhsCSE.cseAdd(newexp, cloc);
                else
                    Add(newcmds, cid(cloc, newexp));
                fi;
            fi;
        fi;
    end,

    procIF := meth(self, c, newcmds, live_out)
        local lo_map, v, orig, then_cmd, else_cmd, then_cp, else_cp;
        c.cond := self.procRHS(c.cond);

        if IsValue(c.cond) then
            Error("This is a constant IF, it should be folded away before reaching copyprop");
            # Following code should work if needed (if you uncomment the error) but it is not optimized
            # self.flush(); # YSV: why is this needed??? I will comment this out.
            if c.cond.v=0 or c.cond.v=false then  Add(newcmds, self.copyProp(c.else_cmd));
            else                                  Add(newcmds, self.copyProp(c.then_cmd));
            fi;
        else
            # Here the idea is that each part of the branch should work with
            # the original csetab (there cannot be crossings) and the final thing
            # should be restored to the original csetab

            then_cp := CopyFields(self, rec(csetab := CopyTab(self.csetab)))
                .prepVarMap(CopyTab(self.varmap), true, false, false);

            else_cp := CopyFields(self, rec(csetab := CopyTab(self.csetab)))
                .prepVarMap(CopyTab(self.varmap), true, false, false);

            then_cmd := then_cp.copyProp(c.then_cmd);
            self.closeVarMap(then_cp, then_cmd);

            else_cmd := else_cp.copyProp(c.else_cmd);
            self.closeVarMap(else_cp, else_cmd);

            Add(newcmds, IF(c.cond, then_cmd, else_cmd));
            #self.flush();
        fi;
    end,

    #if we do Scalar Replacement on the fly, reinject the final writes as assigns (for now)
    finalizeScalarReplacement := meth(self, newcmds)
        local entry;
            if IsBound(self.lhsCSE.csetab.nth) then
                # MRT: sort by order of variable being assigned TO rather
                # than variable assigned FROM. This reduces cache misses by
                # exploiting possible locality in the output array.
                #
                # eg: Y[0]  := a1   ===>  Y[0]   := a1
                #     Y[32] := a2         Y[1]   := a3
                #     Y[1]  := a3         Y[32]  := a2

                for entry in _sortByIdx(self.lhsCSE.csetab.nth) do
                    if Length(entry.args)=2 then
                        Add(newcmds,
                            assign(nth(entry.args[1],entry.args[2]),
                                entry.loc));
                    fi;
                od;
            fi;
            if IsBound(self.lhsCSE.csetab.deref) then
                for entry in self.lhsCSE.csetab.deref do
                    if Length(entry.args)=1 then
                        Add(newcmds,
                            assign(deref(entry.args[1]),
                        entry.loc));
                    fi;
                od;
            fi;
    end,

    # create the reverse mapping, to map last assignment to live_out variable
    #  to the actual variable intead of an SSA related substitute
    substLiveOut := meth(self, newcmds, live_out)
        local varmap, lo_map, active, v;
        varmap := self.varmap;
        lo_map := tab();
        active := Set(Collect(newcmds, var));  # explain this
        for v in live_out do
           # NOTE: The commented out lines with #X removed extra copies
           #  associated with live_out variables. But it turns out that this
           #  does not work if a live_out variable is never actually remapped,
           #  but has an entry in varmap due to copy propagation.
           #  For example (a = t11; b = a) if a is live_out will be compiled to
           #  (b=t11) b CopyPropagate and converted to (b=a) by function below.
           #   which is incorrect

           #X if not IsVar(varmap.(v.id)) or not (varmap.(v.id) in active) then
                Add(newcmds, assign(v, varmap.(v.id)));
           #X else
           #X     lo_map.(varmap.(v.id).id) := v;
           #X fi;
            Unbind(varmap.(v.id));   # prevent duplication of cleanup code due to closeVarMap
        od;
        newcmds := SubstVars(newcmds, lo_map);
        return newcmds;
    end,

    copyProp := meth(self, code)
        local c, cid, cmds, newcmds, live_out;
        cmds := When(ObjId(code)=chain, code.cmds, [code]);
        newcmds := [];
        live_out := Set([]);
        for c in cmds do
            cid := ObjId(c);
            if   IsAssign(c)     then  self.procAssign(c, newcmds, live_out);
            elif (cid = IF)      then  self.procIF(c, newcmds, live_out);
            elif IsExpCommand(c) then  Add(newcmds, self.procRHS(c));
            elif (cid = skip)    then  ;  # do nothing

            # if .sideeffect is bound, the command is assumed have side effects
            # so one has to enable array scalarization inside all its children
            elif IsRec(c) and IsBound(c.sideeffect) and c.sideeffect then
               Add(newcmds, map_children_safe(c, x -> self.procRHS(x)));

            # the command is a container for other functions, copyprop all children
            else Add(newcmds, map_children_safe(c,
			x -> Cond(IsCommand(x), self.copyProp(x), self.procRHS(x))));
            fi;
        od;

        When(self.doScalarReplacement, self.finalizeScalarReplacement(newcmds));
        newcmds := self.substLiveOut(newcmds, live_out);
        # self.flush();
        return chain(newcmds);
    end
));


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

Class(ScriptGenAVX, ScriptGenBase, rec(

	_arch := () -> "AVX",

	_init := meth(self)
		self._setTransform(SGKEY_FFT);
		self._setType(SGKEY_SPCX);
		self._setSize(64);
		self._setFilename("");
	end,

	_validTransforms := meth(arg)
		return [SGKEY_FFT, SGKEY_IFFT, SGKEY_FFT_2D, SGKEY_IFFT_2D];
	end,

	_validSizes := meth(arg)
		local self, szs, xform, type;
		self  := arg[1];
		szs   := [];
		xform := self.getSettingsValue(SGKEY_TRANSFORM);
		type  := self.getSettingsValue(SGKEY_DATATYPE);

		if xform in [SGKEY_FFT, SGKEY_IFFT] then
			if type = SGKEY_DPCX then
				szs := Filtered([1..1024], i->ForAll(Factors(i), j->j<=16) and IsInt(i/16));
				Append(szs, List([11..16], i->2^i));
			else
				szs := Filtered([1..1024], i->ForAll(Factors(i), j->j<=16) and IsInt(i/64));
				Append(szs, List([11..16], i->2^i));
			fi;
		elif xform in [SGKEY_FFT_2D, SGKEY_IFFT_2D] then
			if type = SGKEY_DPCX then
				szs := List([4..8], i -> [2^i, 2^i]);
			else
				szs := List([6..9], i -> [2^i, 2^i]);
			fi;
		elif xform = SGKEY_WHT then
			if type = SGKEY_DPCX then
				szs := List([4..8], i->2^i);
			else
				szs := List([6..8], i->2^i);
			fi;
		fi;

		return szs;
	end,

	_validTypes := meth(arg)
		return [SGKEY_SPCX, SGKEY_DPCX];
	end,

	_genScript := meth(self, runType)
		local scrstr, xform, type, szs, tempvar, tempopts, optrec, optstr, isa, funcstr,
			use_cx, vec2;

		tempvar  := "tdp";
		xform	 := self.getSettingsValue(SGKEY_TRANSFORM);
		type	 := self.getSettingsValue(SGKEY_DATATYPE);
		szs      := self.getSettingsValue(SGKEY_SIZE);
		scrstr	 := "";

		tempopts := tempvar::"Opts";

		if (type = SGKEY_SPCX) then
			isa    := "AVX_8x32f";
			vec2   := 64;
		else
			isa    := "AVX_4x64f";
			vec2   := 16;
		fi;

		use_cx := ForAny(szs, s -> not IsInt(s / vec2));

		if xform in [SGKEY_FFT, SGKEY_IFFT] then
			optrec := rec(
				globalUnrolling    := 128,
				tsplRader          := false,
				tsplBluestein      := false,
				tsplPFA            := false,
				oddSizes           := false,
				interleavedComplex := true,
				cplxVect           := use_cx,
				realVect           := not use_cx,
			);
			if use_cx then
				optrec.RDFT  := false;
            	optrec.URDFT := true;
				if (type = SGKEY_SPCX) then
					optrec.CT      := false;
					optrec.PD      := false;
					optrec.svct    := true;
					optrec.flipIxA := true;
				elif (type = SGKEY_DPCX) then
			        optrec.svct   := false;
					optrec.splitL := true;
				fi;
			fi;
			if xform = SGKEY_IFFT then
				optrec.transInverse := true;
			fi;
			funcstr := "doSimdDft("::String(szs)::", "::isa::", "::tempopts::")";
		elif xform in [SGKEY_FFT_2D, SGKEY_IFFT_2D] then
			optrec := rec(
				interleavedComplex := true,
                oddSizes := false,
				svct := true,
				splitL := false,
				pushTag := true,
				flipIxA := false,
				stdTTensor := true,
				tsplPFA := false
			);
			if xform = SGKEY_IFFT_2D then
				optrec.transInverse := true;
			fi;
			funcstr := "doSimdMddft("::String(szs)::", "::isa::", "::tempopts::")";
		elif xform = SGKEY_WHT then
			optrec := rec(
				verify   := true,
				oddSizes := false,
				svct     := true
			);
			funcstr := "doSimdWht("::String(szs)::", "::isa::", "::tempopts::")";
		else
			return "";
		fi;

		optrec.faultTolerant := true;

		if IsBound(self._settings.(SGKEY_FUNCNAME)) then
			optrec.functionNameRoot := self._settings.(SGKEY_FUNCNAME);
		fi;

		optstr := StringPrint(optrec);
		scrstr := tempopts::" := "::optstr::";;\n";
		Append(scrstr, tempvar::" := "::funcstr::";;\n");
		Append(scrstr, tempvar::"."::runType::"();\n");

		return scrstr;
	end,

	_localSupport := function()
		return (IsBound(LocalConfig) and IsBound(LocalConfig.cpuinfo) and IsBound(LocalConfig.cpuinfo.SIMD().hasAVX) and
			LocalConfig.cpuinfo.SIMD().hasAVX());
	end,

)); # Class ScriptGenAVX

SetScriptGenConstructor(ScriptGenAVX);


text:
local reps,tori,labels,I,II,ords,blocks;

reps:=[[[1,1],[2,1]],
       [[2,1],[5,1]],
       [[4,1],[6,1]],
       [[1,1]],
       [[2,1]]];

I:=[[],[1],[],[],[]];

#
# The parabolic where the representative is distinguished
#
II:=[[1,2],
     [1,2],
     [],
     [1],
     [2]];

labels:=["G_2","G_2(a_1)","(\\tilde A_1)^{(3)}","\\tilde A_1","A_1"];

tori:=[
    [2,2],
    [0,2],
    [],
    [2,-3],
    [-1,2]
    ];

blocks:=[[[9, 1], [5, 1]], [[3, 4], [2, 1]], [[3, 3], [2, 2], [1, 1]], [[3, 3], [1, 5]], [[3, 1], [2, 4], [1, 3]]];

#ords:=[];
ords:=[9,3,3,3,3];

return [reps,tori,labels,I,II,ords,blocks];