text:
################################################################################
##
#W DwG_rank.gi              GroupTheoretical Package
##
#W Paul Bruillard, Cesar Galindo, Siu-Hung Ng, Julia Plavnik, Eric Rowell,
#W Zhenghan Wang
##
## Installation file for DwG_rank functions of the GroupTheoretical Package
##
#Y Copyright (C) 2016, Battelle Memorial Institute
##
################################################################################

###############################################################################
##
#F DwG_rank(<group>,<bool>). . . . . . . .compute rank of D^w(G)
##
InstallGlobalFunction(DwG_rank, function(G,display)
  local A,R,C,CH,H,classes,cs,fs,CC,Counter,f,DwG_rank,x,Zx,ZxReps,alpha_regular,g,is_alpha_regular,Zg_in_Zx,h,alpha_g_h,alpha_h_g,DwG_ranks,last_rank,H3G;

  # compute H_3(G,Z)
  H3G:=GroupHomology(G,3);;
  A:=CyclicGroup(Lcm(H3G));;
  #A:=CyclicGroup(Order(G));
  A:=TrivialGModuleAsGOuterGroup(G,A);;
  R:=ResolutionFiniteGroup(G,4);;
  C:=HomToGModule(R,A);;
  CH:=CohomologyModule(C,3);;
  H:=ActedGroup(CH);
  classes:=Elements(ActedGroup(CH));;
  Print(Length(classes));;
  Print("\n Cohomology group:");;
  Print(StructureDescription(H));;
  Print("\n");

  # get the 3-cocycle rerpresenting the second cohomology classes
  cs:=List(classes,x->CH!.representativeCocycle(x));

  # get a cocycle f:GxGxG->A corresponding to this cohomology class
  fs:=List(cs,x->Mapping(x));

  # get group conjugacy class represenatives
  CC:=List(ConjugacyClasses(G),Representative);;

  # for each cocycle f and each conjugacy class x compute alpha the 2-cocycle on
  # the centralizer
  Counter:=0;
  DwG_ranks:=[];
  for f in fs
  do
    Counter:=Counter+1;
    DwG_rank:=0;
    for x in CC
    do
      # compute the centralizer of X
      Zx:=Centralizer(G,x);;
      ZxReps:=List(ConjugacyClasses(Zx),Representative); # compute the alpha-regular elements of Zx
      alpha_regular:=[];;
      for g in ZxReps
      do
        is_alpha_regular:=true;;
        Zg_in_Zx:=Centralizer(Zx,g);
        for h in Zg_in_Zx
        do
          alpha_g_h:=f(g,h,x)*f(x,g,h)*(f(g,x,h)^-1);;
          alpha_h_g:=f(h,g,x)*f(x,h,g)*(f(h,x,g)^-1);;
          if not(alpha_g_h=alpha_h_g) then
            is_alpha_regular:=false;;
            break;;
          fi;
        od;
        if is_alpha_regular=true then
          Append(alpha_regular,[g]);
          DwG_rank:=DwG_rank+1;
        fi;
      od;
    od;
    if display then
      Print("\nThe cohomology class is: ");
      Print(classes[Counter]);
      Print("\nRank is: ");
      Print(DwG_rank);
    else
      Print(".");
    fi;
    Append(DwG_ranks,[[classes[Counter],DwG_rank]]);
    if Counter=1 then
      last_rank:=DwG_rank;
    fi;
  od;
  Print("\n");

  return DwG_ranks;;
end);

#E DwG_rank.gi . . . . . . . . . . . . . . . . . . . . . . . . . . . . ends here

max_stars_repo_name: pnnl/GroupTheoretical

max_issues_repo_name: pnnl/GroupTheoretical

max_forks_repo_name: pnnl/GroupTheoretical


text:
--
-- The Java Lexer
--
%Options fp=CncLexer
%options single-productions
%options package=CnCParser
%options template=LexerTemplateF.gi
%options filter=cncKWLexer.gi

%Notice
/.
//
// This file is part of the CNC-C implementation and
// distributed under the Modified BSD License.
// See LICENSE for details.
//
// I AM A GENERATED FILE. PLEASE DO NOT CHANGE ME!!!
//
./
%End

%Define
    --
    -- Definition of macro used in the included file LexerBasicMapB.g
    --
    $kw_lexer_class /.$cncKWLexer./
    $_IDENTIFIER    /.$_T_NAME./

%End

%Include
    LexerBasicMapF.gi
%End

%Export

    T_NAME
    T_NUMBER
    T_QUOTEDVAL

    SEMICOLON
    RIGHT_ARROW
    LEFT_ARROW
    COLON_COLON
    COMMA
    LESS_THAN
    GREATER_THAN
    COLON
    LEFT_PARENTHESIS
    RIGHT_PARENTHESIS
    LEFT_BRACKET
    RIGHT_BRACKET
    LEFT_BRACE
    RIGHT_BRACE
    EQUAL
    AMPERSAND
    STAR
    DOT
    DOT_DOT
    PLUS
    MINUS
    SLASH
    SHARP
    ATSIGN
	--UNSIGNED
    --STRUCT

%End

%Terminals
    CtlCharNotWS

    LF   CR   HT   FF

    a    b    c    d    e    f    g    h    i    j    k    l    m
    n    o    p    q    r    s    t    u    v    w    x    y    z
    _

    A    B    C    D    E    F    G    H    I    J    K    L    M
    N    O    P    Q    R    S    T    U    V    W    X    Y    Z

    0    1    2    3    4    5    6    7    8    9

    AfterASCII   ::= '\u0080..\ufffe'
    Space        ::= ' '
    LF           ::= NewLine
    CR           ::= Return
    HT           ::= HorizontalTab
    FF           ::= FormFeed
    DoubleQuote  ::= '"'
    SingleQuote  ::= "'"
    Percent      ::= '%'
    VerticalBar  ::= '|'
    Exclamation  ::= '!'
    AtSign       ::= '@'
    BackQuote    ::= '`'
    Tilde        ::= '~'
    Sharp        ::= '#'
    DollarSign   ::= '$'
    Ampersand    ::= '&'
    Caret        ::= '^'
    Colon        ::= ':'
    SemiColon    ::= ';'
    BackSlash    ::= '\'
    LeftBrace    ::= '{'
    RightBrace   ::= '}'
    LeftBracket  ::= '['
    RightBracket ::= ']'
    QuestionMark ::= '?'
    Comma        ::= ','
    Dot          ::= '.'
    LessThan     ::= '<'
    GreaterThan  ::= '>'
    Plus         ::= '+'
    Minus        ::= '-'
    Slash        ::= '/'
    Star         ::= '*'
    LeftParen    ::= '('
    RightParen   ::= ')'
    Equal        ::= '='
    -- Unsigned	 ::= "unsigned"
    -- Struct		 ::= "struct"

%End

%Start
    Token
%End

%Rules

    ---------------------  Rules for Scanned Tokens --------------------------------
    -- The lexer creates an array list of tokens which is defined in the PrsStream class.
    -- A token has three attributes: a start offset, an end offset and a kind.
    --
    -- Only rules that produce complete tokens have actions to create token objects.
    -- When making a token, calls to the methods, $getToken(1) and $getRightSpan(),
    -- provide the offsets (i.e. the span) of a rule's right hand side (rhs) and thus of the token.
    -- For a rule of the form A ::= A1 A2 ... An, the start offset of the rhs of A is given by
    -- $getToken(1) or by $getLeftSpan() and the end offset by $getRightSpan().
    --
    -- Regarding rules for parsing in general, note that for a rhs symbol Ai, the
    -- method $getToken(i) returns the location of the leftmost character derived from Ai.
    -- The method $getLeftSpan(i) returns the same location unless Ai produces %Empty in which case
    -- it returns the location of the last character derived before reducing Ai to %Empty.
    -- The method $getRightSpan(i) returns the location of the rightmost character derived from Ai
    -- unless Ai produces %Empty in which case it returns the location of the last character
    -- derived before reducing Ai to %Empty.
    --------------------------------------------------------------------------------
    Token ::= ';'
        /.
                    makeToken($_SEMICOLON);
        ./

    Token ::= '-' '>'
        /.
                    makeToken($_RIGHT_ARROW);
        ./

    Token ::= '<' '-'
        /.
                    makeToken($_LEFT_ARROW);
        ./

    Token ::= ':' ':'
        /.
                    makeToken($_COLON_COLON);
        ./

    Token ::= ','
        /.
                    makeToken($_COMMA);
        ./

    Token ::= '<'
        /.
                    makeToken($_LESS_THAN);
        ./

    Token ::= '>'
        /.
                    makeToken($_GREATER_THAN);
        ./

    Token ::= ':'
        /.
                    makeToken($_COLON);
        ./

    Token ::= '('
        /.
                    makeToken($_LEFT_PARENTHESIS);
        ./

    Token ::= ')'
        /.
                    makeToken($_RIGHT_PARENTHESIS);
        ./

    Token ::= '['
        /.
                    makeToken($_LEFT_BRACKET);
        ./

    Token ::= ']'
        /.
                    makeToken($_RIGHT_BRACKET);
        ./

	Token ::= '{'
        /.
                    makeToken($_LEFT_BRACE);
        ./

    Token ::= '}'
        /.
                    makeToken($_RIGHT_BRACE);
        ./

    Token ::= '='
        /.
                    makeToken($_EQUAL);
        ./

	Token ::= '@'
        /.
                    makeToken($_ATSIGN);
        ./

    Token ::= Identifier
        /.
                    checkForKeyWord();
        ./

    Token ::= '&'
        /.
                    makeToken($_AMPERSAND);
        ./

    Token ::= '*'
        /.
                    makeToken($_STAR);
        ./

    Token ::= '/'
        /.
                    makeToken($_SLASH);
        ./

	Token ::= '+'
        /.
                    makeToken($_PLUS);
        ./

    Token ::= '-'
        /.
                    makeToken($_MINUS);
        ./

    Token ::= '.'
        /.
                    makeToken($_DOT);
        ./

	Token ::= '.' '.'
        /.
                    makeToken($_DOT_DOT);
        ./

   	Token ::= '#'
        /.
                    makeToken($_SHARP);
        ./
    --Token ::= 'u' 'n' 's' 'i' 'g' 'n' 'e' 'd'
      --  /.
       --             makeToken($_UNSIGNED);
       -- ./
    --Token ::= 's' 't' 'r' 'u' 'c' 't'
      --  /.
        --            makeToken($_STRUCT);
        --./

    Token ::= '"' SLBody '"'
        /.
                    makeToken($_T_QUOTEDVAL);
        ./

    Token ::= '"' SLBody Eol -- TODO: remove the Eol from the string
        /.
                    makeToken($_T_QUOTEDVAL);
        ./

    Token ::= Integer
        /.
                    makeToken($_T_NUMBER);
        ./
    Token ::= WS -- White Space is scanned but not added to output vector
        /.
                    skipToken();
        ./

    Token ::= SLC
        /.
                    skipToken();
        ./

    SLBody ::= %Empty
             | SLBody NotDQ

    SLC ::= '/' '/'
          | SLC NotEol

    NotEol -> Letter
            | Digit
            | Special
            | Space
            | HT
            | FF
            | CtlCharNotWS

    Special -> '+' | '-' | '*' | '/' | '(' | ')' | '"' | '!' | '@' | '`' | '~' |
               '%' | '&' | '^' | ':' | ';' | "'" | '\' | '|' | '{' | '}' |
               '[' | ']' | '?' | ',' | '.' | '<' | '>' | '=' | '#'

    SpecialNotDQ -> '+' | '-' | '*' | '/' | '(' | ')' | '!' | '@' | '`' | '~' |
                    '%' | '&' | '^' | ':' | ';' | "'" | '\' | '|' | '{' | '}' |
                    '[' | ']' | '?' | ',' | '.' | '<' | '>' | '=' | '#'

    Integer -> Digit
             | Integer Digit

    WSChar -> Space
            | Eol
            | HT
            | FF

    Eol -> LF
         | CR

    Letter -> LowerCaseLetter
            | UpperCaseLetter
            | _
            | '$'
            | '\u0080..\ufffe'

    LowerCaseLetter -> a | b | c | d | e | f | g | h | i | j | k | l | m |
                       n | o | p | q | r | s | t | u | v | w | x | y | z

    UpperCaseLetter -> A | B | C | D | E | F | G | H | I | J | K | L | M |
                       N | O | P | Q | R | S | T | U | V | W | X | Y | Z

    Digit -> 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9

    WS -> WSChar
        | WS WSChar

    Identifier -> Letter
                | Identifier Letter
                | Identifier Digit

    NotDQ -> Letter
           | Digit
           | SpecialNotDQ
           | Space
           | HT
           | FF
%End

ext: gi

max_stars_repo_name: pelmers/cnc-ocr

max_issues_repo_name: pelmers/cnc-ocr

max_forks_repo_name: pelmers/cnc-ocr


text:
Composition := function(f, g)
    return x -> f(g(x));
end;

h := Composition(x -> x+1, x -> x*x);
h(5);
# 26

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

Declare(compareMatrices);
Declare(applyCellRules);
Declare(doCellWht);

testfunc := function(m1)
   m1.isRight := true;
   return(m1);
end;

subtractDFTMatrices := function(m1, m2)
   local mdiff, verified;

    m1 := Cond(
            ObjId(m1) = ListClass,
            m1,
            IsBound(m1.isSPL) and m1.isSPL,
            MatSPL(m1),
            Error("Don't know what m1 is"));
    m2 := Cond(
            ObjId(m2) = ListClass,
            m2,
            IsBound(m2.isSPL) and m2.isSPL,
            MatSPL(m2),
            Error("Don't know what m2 is"));

return Maximum(Flat(MapMat(m1-m2, i->AbsComplex(ComplexAny(i)))));

end;

compareMatrices := function(m1, m2)
   local mdiff, verified;

    m1 := Cond(
            ObjId(m1) = ListClass,
            m1,
            IsBound(m1.isSPL) and m1.isSPL,
            MatSPL(m1),
            Error("Don't know what m1 is"));
    m2 := Cond(
            ObjId(m2) = ListClass,
            m2,
            IsBound(m2.isSPL) and m2.isSPL,
            MatSPL(m2),
            Error("Don't know what m2 is"));

    mdiff := m1-m2;

    #return(Cond(Minimum(Flat(mdiff)) = Maximum(Flat(mdiff)), true, false));
    verified := InfinityNormMat(mdiff);
    When(verified <= 0.001, Print(gap.colors.Green("Ok."), "\n"), gap.colors.Red("Matrices are DIFFERENT!\n"));
    return(verified);
end;

# n := 4; spus := 2; pkSize := 1;
doCellWht := function(n, spus, pkSize, opts)
    local t, r, sums, sumsorig, c, m, mdiff, measured, verified, compare;

    opts.spus := spus;

    t       := WHT(n).withTags([ ParCell(spus, pkSize) ]);
    r       := RandomRuleTree(t, opts);
    sums    := SumsRuleTreeOpts(r, opts);
    #sumsorig:= Copy(sums);
    #sums    := applyCellRules(sums, opts);
    #Print("Comparing sums after application of rules...");
    #compareMatrices(sums, sumsorig);
    if compareMatrices(sums, t) = false then
      return([r, sums, ""]);
    fi;

    c := CodeSumsOpts(sums, opts);
    measured := CMeasure(c, opts);
    Print(measured, " [ cycles ]\n");
    m := CMatrix(c, opts);
    Print("Comparing via CMatrix verification...");
    compareMatrices(m, t);

    return([r, sums, m]);

end;


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

ImportAll(paradigms.smp);
ImportAll(paradigms.vector);

Class(FFTXPOWER9Opts, FFTXOpts, rec(
    tags := [],
    operations := rec(Print := s -> Print("<FFTX POWER9 options record>")),
));

power9Opts := function(arg)
    local opts, optrec, vsxopts;

    optrec := rec(dataType := T_Real(64), globalUnrolling := 32);
    vsxopts := rec(svct := true, splitL := false, oddSizes := false, stdTTensor := true, tsplPFA := false, realVect := false, cplxVect := true);

    opts := CopyFields(FFTXOpts, SIMDGlobals.getOpts(POWER9_2xf, vsxopts));
    opts.breakdownRules.TFCall := FFTXOpts.breakdownRules.TFCall;

    opts.globalUnrolling := optrec.globalUnrolling;

    # FFTX specific breakdown rules
    opts.breakdownRules.Circulant := [Circulant_PRDFT_FDataNT];
    opts.breakdownRules.PRDFT := List([PRDFT1_Base1, PRDFT1_Base2, PRDFT1_CT, PRDFT1_PF, PRDFT_PD, PRDFT_Rader], _noT);
    opts.breakdownRules.IPRDFT := List([ IPRDFT1_Base1, IPRDFT1_Base2, IPRDFT1_CT, IPRDFT_PD, IPRDFT_Rader ], _noT);
    opts.breakdownRules.PRDFT3 := List([ PRDFT3_Base1, PRDFT3_Base2, PRDFT3_CT ], _noT);

    return opts;
end;

Declare(ParseOptsPOWER);

Class(FFTXPOWER9DefaultConf, rec(
    __call__ := self >> self,
    getOpts := (self, t) >> ParseOptsPOWER(self, t),
    operations := rec(Print := s -> Print("<FFTX POWER9 Default Configuration>")),
    useOMP := false,
    useSIMD := true
));

power9Conf := rec(
    defaultName := "defaultPOWER9Conf",
    defaultOpts := (arg) >> FFTXPOWER9DefaultConf,
    useOMP := false,
    useSIMD := true,
    confHandler := power9Opts
);

fftx.FFTXGlobals.registerConf(power9Conf);

getTargetOS := function()
    local tgt;

    if LocalConfig.osinfo.isWindows() then
        tgt := "win-x64-cuda";
    elif LocalConfig.osinfo.isLinux() then
        tgt := "linux-cuda";
    elif LocalConfig.osinfo.isDarwin() then
        tgt := "linux-cuda";    ## may work
    fi;
    return tgt;
end;

#-----------------------------------------
# OpenMP + VSX

Class(FFTXPOWER9OMPOpts, FFTXOpts, rec(
    tags := [],
    operations := rec(Print := s -> Print("<FFTX POWER9 OpenMP options record>")),
));

power9OMPOpts := function(arg)
    local opts, optrec, vsxopts, smpopts, tid;

    optrec := rec(dataType := T_Real(64), globalUnrolling := 32);
    vsxopts := rec(svct := true, splitL := false, oddSizes := false, stdTTensor := true, tsplPFA := false);
    smpopts := rec(numproc := LocalConfig.cpuinfo.cores, api := "OpenMP");

    opts := CopyFields(FFTXOpts, SIMDGlobals.getOpts(POWER9_2xf, vsxopts));
    opts.breakdownRules.TFCall := FFTXOpts.breakdownRules.TFCall;

    #-- OpenMP opts merging--
    opts.unparser := When(IsBound(smpopts.OmpMode) and smpopts.OmpMode = "for",
                            OpenMP_POWERUnparser_ParFor,
                            OpenMP_POWERUnparser);
    opts.codegen := spiral.libgen.VecRecCodegen;

    opts.breakdownRules.GT := [ GT_Base, GT_NthLoop,
        CopyFields(GT_Par, rec(parEntireLoop := false, splitLoop := true)), GT_Par_odd,
        GT_Vec_AxI, GT_Vec_IxA, GT_Vec_IxA_Push, GT_Vec_IxA_L, GT_Vec_L_IxA
    ];
    opts.breakdownRules.TTensorI := Concat([
    CopyFields(TTensorI_toGT, rec(
        applicable := (self, t) >> t.hasTags() and ObjId(t.getTags()[1])=AParSMP ))],
        opts.breakdownRules.TTensorI);

    opts.breakdownRules.TTensorInd :=
        Concat([dsA_base_smp, dsA_smp, L_dsA_L_base_smp, L_dsA_L_smp],
            opts.breakdownRules.TTensorInd);

    tid := When(smpopts.api = "OpenMP", threadId(), var("tid", TInt));
    opts.tags := Concat([ AParSMP(smpopts.numproc, tid)  ], opts.tags);
    #----------
    opts.globalUnrolling := optrec.globalUnrolling;

    # FFTX specific breakdown rules
    opts.breakdownRules.Circulant := [Circulant_PRDFT_FDataNT];
    opts.breakdownRules.PRDFT := List([PRDFT1_Base1, PRDFT1_Base2, PRDFT1_CT, PRDFT1_PF, PRDFT_PD, PRDFT_Rader], _noT);
    opts.breakdownRules.IPRDFT := List([ IPRDFT1_Base1, IPRDFT1_Base2, IPRDFT1_CT, IPRDFT_PD, IPRDFT_Rader ], _noT);
    opts.breakdownRules.PRDFT3 := List([ PRDFT3_Base1, PRDFT3_Base2, PRDFT3_CT ], _noT);

    return opts;
end;

Declare(ParseOptsPOWER);

Class(FFTXPOWER9OMPConf, rec(
    __call__ := self >> self,
    getOpts := (self, t) >> ParseOptsPOWER(self, t),
    operations := rec(Print := s -> Print("<FFTX POWER9 OpenMP Configuration>")),
    useOMP := true,
    useSIMD := true
));

power9OMPConf := rec(
    defaultName := "defaultPOWER9OMPConf",
    defaultOpts := (arg) >> FFTXPOWER9OMPConf,
    useOMP := true,
    useSIMD := true,
    confHandler := power9OMPOpts
);

fftx.FFTXGlobals.registerConf(power9OMPConf);

getTargetOS := function()
    local tgt;

    if LocalConfig.osinfo.isWindows() then
        tgt := "win-x64-cuda";
    elif LocalConfig.osinfo.isLinux() then
        tgt := "linux-cuda";
    elif LocalConfig.osinfo.isDarwin() then
        tgt := "linux-cuda";    ## may work
    fi;
    return tgt;
end;

# this is a first experimental opts-deriving logic. This needs to be done extensible and properly
ParseOptsPOWER := function(conf, t)
    local tt, _tt, _conf, _opts;

    if conf.useOMP then
        _opts := power9OMPOpts();
    else
        _opts := power9Opts();
    fi;
    return _opts;

    Error("Don't know how to derive opts!\n");
end;

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:
# Copyright (c) 2018-2021, Carnegie Mellon University
# See LICENSE for details

Class(RulesHfunc, RuleSet, rec(__avoid__ := [])); #Diag]));

# IsInt(range(e)) makes sure that we don't apply this to functions with range of TInt,
# which are special. Otherwise
# the following breaks eventually
# RC(Diag(fPrecompute(fCompose(dOmega(6, 1), diagTensor(fCompose(dLin(V(2), 1, 0, TInt),
#                       fBase(i601)), dLin(3, 1, 0, TInt)), H(3, 1, 0, 1))))),
# because this rule matches
base := pat -> pat.cond(e->domain(e)=1 and IsInt(range(e)));

unmemo := exp -> When(IsVar(exp) and IsBound(exp.mapping), exp.mapping, exp);
SetRange := (exp,r) >> When(IsVar(exp), exp.setRange(r), exp);

RewriteRules(RulesHfunc, rec(
 # fId -> H
 #Hid := Rule(fId, e -> H(e.params[1], e.params[1], 0, 1)),

 # fTensor -> H o f
 Tensor_toH1 := ARule(fTensor, [ @F, base(@(1)) ],
     (e,cx) -> [ fCompose(H(range(@(1).val) * range(@F.val),
		            range(@F.val),
		            nomemo(cx, "b", @(1).val.lambda().at(0)),
			    range(@(1).val)),
	                  @F.val) ]),
 # fTensor -> H o f
 Tensor_toH2 := ARule(fTensor, [ base(@(1)), @F ],
     (e,cx) -> [ fCompose(H(range(@(1).val)*range(@F.val),
		            range(@F.val),
		            nomemo(cx, "b", range(@F.val) * @(1).val.lambda().at(0)),
			    1),
	                  @F.val) ]),
 # gammaTensor -> HZ o f
 GammaTensor_toHZ1 := ARule(gammaTensor, [ @F, base(@(1)) ],
     (e, cx) -> let(F := @F.val, b := @(1).val, j := b.lambda().at(0),
	 [ fCompose(HZ(range(b)*range(F), range(F), memo(cx,"g",range(F)*j), range(b)), F) ])),

 # gammaTensor -> HZ o f, same RHS as in rule above
 GammaTensor_toHZ2 := ARule(gammaTensor, [ base(@(1)), @F ], #~.GammaTensor_toHZ1[3]),
     (e, cx) -> let(F := @F.val, b := @(1).val, j := b.lambda().at(0),
	 [ fCompose(HZ(range(b)*range(F), range(F), memo(cx,"g",range(F)*j), range(b)), F) ])),

 # H o H -> H
 H_H := ARule(fCompose, [ [H, @N, @n, @b, @s], [H, @n, @m, @bb, @ss] ],
     (e, cx) -> [ H(@N.val, @m.val,
	            memo(cx, "b", @b.val + @s.val * @bb.val), # new base
		    @s.val * @ss.val) ]),                     # new stride
 ## HZ o (H | HZ) -> HZ
 HZ_H := ARule(fCompose, [ [HZ, @N, @n, @b, @s], [@(1,[H,HZ]), @n, @m, @bb, @ss] ],
     (e, cx) -> [ HZ(@N.val, @m.val,
	             memo(cx, "g", @b.val + @s.val * @bb.val), # new base
		     @s.val * @ss.val) ]),                     # new stride
 # RM o (H | HZ) -> RM  ("rho" in the PLDI paper)
 RM_H := ARule(fCompose, [ [RM, @N, @n, @phi, @g], [@(1,[HZ,H]), @n, @m, @b, @s] ],
     (e, cx) -> [ RM(@N.val, @m.val,
	             SetRange(memo(cx, "f", powmod(@phi.val,@g.val,unmemo(@b.val),@N.val)), @N.val),
		     powmod(1, @g.val, @s.val, @N.val).eval()) ]),

# Following rule is for use in expressions like Gath(fBase(i))
 H_fBase := Rule(@(1,fBase),e ->H(@(1).val.params[1],1,@(1).val.params[2],1)),

#Handle (H tensor I) o H
# H_tensorI := ARule(fCompose,[[fTensor,[H,@(1),@(2),@(3),@(4)],[fId,@(5)]],[H,@.cond(e->e=@(2).val*@(5).val),1,@(6),@]],
#e->[H(@(1).val*@(5).val,1,(@(3).val+@(4).val*idiv(@(6).val,@(5).val))*@(5).val+imod(@(6).val,@(5).val),1)]),

));


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

# -------------------------------------------------------------------------
# Scalar bit-level instructions
# 16x1i: 16-way bit register
# -------------------------------------------------------------------------

Class(sklr_bcast_16x1i,  VecExp_16.binary());
Class(sklr_loadu_16x1i,  VecExp_16.ternary());
Class(sklr_storeu_16x1i, VecExpCommand.quad(), rec(
    isStoreop := True,
    __call__ := (self, loc, offs, exp, p) >> WithBases(self,
       rec(operations := CmdOps,
       loc := toAssignTarget(loc),
       offs := toExpArg(offs),
       exp := toExpArg(exp),
       p := p)),

    rChildren := self >> [self.loc, self.offs, self.exp, self.p],
    rSetChild := rSetChildFields("loc", "offs", "exp", "p"),

    print := (self,i,is) >> Print(self.__name__, "(", self.loc.print(), ", ", self.offs.print(), ", ", self.exp.print(), ", ", self.p, ")"),

    # in case of explicit type cast (YSV modification), we don't need getNoScalar,
    # and below returns [], compiler understands not to mess with typecasts
    getNoScalar := self >> When(IsBound(self.noscalar) and IsBound(self.args[self.noScalar].loc),
    self.args[self.noScalar].loc, [])
));

# -------------------------------------------------------------------------
# Scalar bit-level instructions
# 32x1i: 32-way bit register
# -------------------------------------------------------------------------

Class(sklr_bcast_32x1i,  VecExp_32.binary());
Class(sklr_loadu_32x1i,  VecExp_32.ternary());
Class(sklr_storeu_32x1i, VecExpCommand.quad(), rec(
    isStoreop := True,
    __call__ := (self, loc, offs, exp, p) >> WithBases(self,
       rec(operations := CmdOps,
       loc := toAssignTarget(loc),
       offs := toExpArg(offs),
       exp := toExpArg(exp),
       p := toExpArg(p))),

    rChildren := self >> [self.loc, self.offs, self.exp, self.p],
    rSetChild := rSetChildFields("loc", "offs", "exp", "p"),

    print := (self,i,is) >> Print(self.__name__, "(", self.loc.print(), ", ", self.offs.print(), ", ", self.exp.print(), ", ", self.p, ")"),

    # in case of explicit type cast (YSV modification), we don't need getNoScalar,
    # and below returns [], compiler understands not to mess with typecasts
    getNoScalar := self >> When(IsBound(self.noscalar) and IsBound(self.args[self.noScalar].loc),
    self.args[self.noScalar].loc, [])
));

# -------------------------------------------------------------------------
# Scalar bit-level instructions
# 64x1i: 64-way bit register
# -------------------------------------------------------------------------

Class(sklr_bcast_64x1i, VecExp_64.binary());
Class(sklr_loadu_64x1i, VecExp_64.ternary());
Class(sklr_storeu_64x1i, VecExpCommand.quad(), rec(
    isStoreop := True,
    __call__ := (self, loc, offs, exp, p) >> WithBases(self,
       rec(operations := CmdOps,
       loc := toAssignTarget(loc),
       offs := toExpArg(offs),
       exp := toExpArg(exp),
       p := toExpArg(p))),

    rChildren := self >> [self.loc, self.offs, self.exp, self.p],
    rSetChild := rSetChildFields("loc", "offs", "exp", "p"),

    print := (self,i,is) >> Print(self.__name__, "(", self.loc.print(), ", ", self.offs.print(), ", ", self.exp.print(), ", ", self.p, ")"),

    # in case of explicit type cast (YSV modification), we don't need getNoScalar,
    # and below returns [], compiler understands not to mess with typecasts
    getNoScalar := self >> When(IsBound(self.noscalar) and IsBound(self.args[self.noScalar].loc),
    self.args[self.noScalar].loc, [])
));


text:
#############################################################################
##
#W    matgrp.gi                                                  Bettina Eick
##
##    This file contains the header functions to handle the 3- and
##    4-dimensional almost crystallographic integral matrix groups.
##

#############################################################################
##
#F AlmostCrystallographicDim3( type, param )
##
InstallGlobalFunction( AlmostCrystallographicDim3,
function( type, param )
    local g, G, info;

    # type is integer or string
    if IsString( type ) then
        g := Int( type );
    elif IsInt( type ) then
        g := type;
    else
        g := false;
    fi;
    if not g in [1..17] then
        Error("type does not define a valid ac-group type");
    fi;

    # check parameters
    if IsBool( param ) then
        param := List( [1..ACDim3Param[g]], x -> Random(Integers) );
        if param[1] = 0 then param[1] := Random(Integers)^2 + 1; fi;
    elif IsInt( param ) then
        param := List( [1..ACDim3Param[g]], x -> param );
    elif Length( param ) <> ACDim3Param[g] then
        Error("parameter should be a list of length ", ACDim3Param[g] );
    fi;

    # get group
    if Length( param ) = 4 then
        G := ACDim3Funcs[g](param[1], param[2], param[3], param[4]);
    elif Length(param) = 2 then
        G := ACDim3Funcs[g](param[1], param[2] );
    elif Length(param) = 1 then
        G := ACDim3Funcs[g](param[1]);
    fi;

    # set some information
    info := rec( dim := 3, type := g, param := param );
    SetAlmostCrystallographicInfo( G, info );
    SetIsAlmostCrystallographic( G, true );
    SetDimensionOfMatrixGroup( G, 4 );
    SetFieldOfMatrixGroup( G, Rationals );
    SetSize( G, infinity );
    return G;
end );

#############################################################################
##
#F AlmostCrystallographicDim4( type, param )
##
InstallGlobalFunction( AlmostCrystallographicDim4,
function( type, param )
    local g, G, info;

    # type is integer or string
    if IsString( type ) then
        g := Position( ACDim4Types, type );
    elif IsInt( type ) then
        g := type;
    else
        g := false;
    fi;
    if not g in [1..95] then
        Error("type does not define a valid ac-group type");
    fi;

    # check parameters
    if IsBool( param ) then
        param := List( [1..ACDim4Param[g]], x -> Random(Integers) );
        if param[1] = 0 then param[1] := Random(Integers)^2 + 1; fi;
    elif IsInt( param ) then
        param := List( [1..ACDim4Param[g]], x -> param );
    elif Length( param ) <> ACDim4Param[g] then
        Error("parameter should be a list of length ", ACDim4Param[g] );
    fi;

    # get group
    if Length(param) = 7 then
        G := ACDim4Funcs[g]( param[1], param[2], param[3], param[4],
                             param[5], param[6], param[7]);
    elif Length(param) = 6 then
        G := ACDim4Funcs[g]( param[1], param[2], param[3], param[4],
                             param[5], param[6]);
    elif Length(param) = 5 then
        G := ACDim4Funcs[g]( param[1], param[2], param[3], param[4],
                             param[5]);
    elif Length(param) = 4 then
        G := ACDim4Funcs[g]( param[1], param[2], param[3], param[4]);
    elif Length(param) = 3 then
        G := ACDim4Funcs[g]( param[1], param[2], param[3]);
    fi;

    # set some information
    info := rec( dim := 4, type := g, param := param );
    SetAlmostCrystallographicInfo( G, info );
    SetIsAlmostCrystallographic( G, true );
    SetDimensionOfMatrixGroup( G, 5 );
    SetFieldOfMatrixGroup( G, Rationals );
    SetSize( G, infinity );
    return G;
end );

#############################################################################
##
#F AlmostCrystallographicGroup( dim, type, param )
##
InstallGlobalFunction( AlmostCrystallographicGroup,
function( dim, type, param )
    if dim = 3 then
        return AlmostCrystallographicDim3( type, param );
    elif dim = 4 then
        return AlmostCrystallographicDim4( type, param );
    else
        Error("dimension must be 3 or 4");
    fi;
end );

#############################################################################
##
#F IsAlmostCrystallographic( <G> )
##
InstallMethod( IsAlmostCrystallographic, "for groups", true,
    [IsGroup ], 0,
function( G )
    if HasAlmostCrystallographicInfo( G ) then
        return true;
    else
        Print("sorry - cannot check this property \n");
        return fail;
    fi;
end );

max_stars_repo_name: alex-konovalov/aclib

max_issues_repo_name: alex-konovalov/aclib

max_forks_repo_name: alex-konovalov/aclib


text:
#############################################################################
##
#W  mindeg.gi            GAP 4 package SingerAlg                Thomas Breuer
##
##  This file contains the implementation for GAP functions related to
##  the numbers <M>m(q,e)</M>.
##

#############################################################################
##
#M  SingerAlgE(  )
#M  SingerAlgE( <q>[, <n>], <z> )
##
InstallMethod( SingerAlgE,
    [ "IsSingerAlgebra" ],
    A -> CallFuncList( SingerAlgE, ParametersOfSingerAlgebra( A ) ) );

InstallMethod( SingerAlgE,
    [ "IsPosInt", "IsPosInt", "IsPosInt" ],
    { q, n, z } -> (q^n-1) / z );

InstallMethod( SingerAlgE,
    [ "IsPosInt", "IsPosInt" ],
    { q, z } -> (q^OrderMod( q, z )-1) / z );

#############################################################################
##
#M  MinimalDegreeOfSingerAlgebraGAP(  )
#M  MinimalDegreeOfSingerAlgebraGAP( <q>, <e> )
#M  MinimalDegreeOfSingerAlgebraGAP( <q>, <n>, <e> )
##
##  If the <Ref Attr="LoewyStructureInfo" Label="for a Singer algebra"/>
##  value of the Singer algebra A is known then use it.
##  If the cheap criteria suffice then take their result.
##  Otherwise call the hard method for q and e;
##  if the dimension of the algebra is small then use the algebra
##  also if just the parameters q and e are given.
##
InstallMethod( MinimalDegreeOfSingerAlgebraGAP,
    [ "IsSingerAlgebra and HasLoewyStructureInfoGAP" ],
    A -> LoewyStructureInfoGAP( A ).m );

InstallMethod( MinimalDegreeOfSingerAlgebraGAP,
    [ "IsSingerAlgebra" ],
    function( A )
    local paras, q, n, e;

    # Note that we cannot assume that the parameter 'n' is equal to
    # the multiplicative order of 'q' modulo 'e'.
    paras:= ParametersOfSingerAlgebra( A );
    q:= paras[1];
    n:= paras[2];
    e:= SingerAlgE( A );
    return MinimalDegreeOfSingerAlgebraGAP( q, OrderModExt( q, e, n ), e );
    end );

InstallMethod( MinimalDegreeOfSingerAlgebraGAP,
    [ "IsPosInt", "IsPosInt" ],
    { q, e } -> MinimalDegreeOfSingerAlgebraGAP( q, OrderMod( q, e ), e ) );

InstallMethod( MinimalDegreeOfSingerAlgebraGAP,
    [ "IsPosInt", "IsPosInt", "IsPosInt" ],
    function( q, n, e )
    local e_str, cache_e, m;

    q:= q mod e;
    e_str:= String( e );  # large integers cannot be record component names

    if not IsBound( MinimalDegreeOfSingerAlgebraCache.( e_str ) ) then
      cache_e:= rec();
      MinimalDegreeOfSingerAlgebraCache.( e_str ):= cache_e;
    else
      cache_e:= MinimalDegreeOfSingerAlgebraCache.( e_str );
      if IsBound( cache_e.( q ) ) then
        return cache_e.( q );
      fi;
    fi;

    m:= MinimalDegreeCheapGAP( q, n, e );
    if m = 0 then
      m:= MinimalDegreeHardGAP( q, n, e );
    fi;

    # Extend the cache.
    cache_e.( q ):= m;

    return m;
    end );

#############################################################################
##
#F  MinimalDegreeCheapGAP( <q>, <n>, <e> )
##
##  Return either <C>0</C> or the number <M>m( q, e )</M>,
##  which is the smallest number of powers of q
##  such that e divides the sum of these powers.
##

##  The return value <C>0</C> means that the cheap criteria from
##  the two papers  and
##  do not suffice to determine <M>m( q, e )</M>.
##

##  We assume that
##  q is smaller than e,
##  the argument n is
##  <C>OrderMod( </C>q<C>, </C>e<C> )</C>.
##  The function does not use cached values of <M>m( q, e )</M>
##  or values known via the database of low dimensional Singer algebras.
##
InstallGlobalFunction( MinimalDegreeCheapGAP, function( q, n, e )
    local e1, e2, m;

    if e = 1 then
      return 1;
    fi;

    Assert( 2, PowerMod( q, n, e ) = 1 and OrderModExt( q, e, n ) = n );

    if q = 1 then
      # Example I.6.1.
      return e;
    elif IsEvenInt( n ) and PowerModInt( q, n/2, e ) = e-1 then
      # Decide the case m = 2, see Lemma I.6.3.
      return 2;
    elif n = 2 then
      # Prop. II.2.7.
      e1:= Gcd( e, q-1 );
      e2:= Gcd( e, q+1 );
      if e2 <= e1 or ( IsEvenInt( e ) and IsEvenInt( (q^2-1)/e ) ) then
        return e1;
      else
        return 2 * e1;
      fi;
    fi;

    # Deal with those cases where e is a prime power,
    # and where we know the value of m(q,e).
    if IsPrimePowerInt( e ) then
      # Let e = p^k for a prime p.
      if e mod 2 = 0 then
        # see Prop. II.2.9, note that here we have e > 4 and q mod e <> q-1.
        if q mod 4 = 1 then
          return Gcd( e, q-1 );
        else
          Assert( 1, ( q+1 ) mod e <> 0 );
          return 4;
        fi;
      else
        # see Prop. II.2.10, note that q \equiv 1 \pmod{p}
        # if and only if \gcd( e, q-1 ) \not= 1 holds.
        m:= Gcd( e, q-1 );
        if m <> 1 then
          return m;
        elif n mod 3 = 0 then
          # Here we know that m(q,e) > 2 and thus p \not= 3.
          # Hence \gcd( e, q^{n/3}-1 ) = 1 holds,
          # because otherwise q^{n/3} would be a p-element in (Z/eZ)^*.
          # Thus Lemma I.6.4 yields m(q,e) \leq 3.
          # (This is a generalization of Cor. II.2.11.)
          return 3;
        elif 2 * n = Phi( e ) then
          # Cor. II.2.14.
          return 3;
        elif e mod 11 = 0 then
          # Prop. II.2.16.
          return 5;
        fi;
      fi;
    fi;

    # Deal with those cases where e is twice an odd prime power,
    # and where we know the value of m(q,e).
    if IsEvenInt( e ) and IsPrimePowerInt( e/2 ) then
      # see Prop. II.2.17, note that the order of q mod e is a power of p
      # if and only if \gcd( e/2, q-1 ) \not= 1 holds.
      m:= Gcd( e, q-1 );
      if m > 2 then
        return m;
      fi;
    fi;

    # We give up.
    return 0;
    end );

#############################################################################
##
#F  VectorIterator( m, n, s )
##
##  iterator for all vectors of length n
##  with entries in { 0, 1, 2, ... m }
##  and coefficient sum s, such that ...
##
##  fill the first n positions in v with smallest partition of s into
##  numbers at most m
##  (we assume that s <= m*n)
##
BindGlobal( "VectorIterator_ResetPrefix", function( v, m, n, s )
    local rest, i, j;

    rest:= s;
    i:= 1;
    while m < rest do
      v[i]:= m;
      rest:= rest - m;
      i:= i + 1;
    od;
    v[i]:= rest;
    for j in [ i+1 .. n ] do
      v[j]:= 0;
    od;

    return v;
    end );

BindGlobal( "VectorIterator", function( m, n, s )
    local next;

    if s > m*n then
      # The iterator is empty.
      next:= false;
    else
      # Initialize the coefficient vector.
      next:= VectorIterator_ResetPrefix( [], m, n, s );
    fi;

    return IteratorByFunctions( rec(
             NextIterator:= function( iter )
               local result, pos, succ;

               result:= iter!.next;
               # Find the first position with a nonzero value
               # such that the value on the right is smaller than m.
               pos:= First( [ 1 .. iter!.n - 1 ],
                            i -> result[i] <> 0 and result[ i+1 ] < iter!.m );
               if pos = fail then
                 succ:= false;
               else
                 succ:= ShallowCopy( result );
                 succ[ pos ]:= result[ pos ] - 1;
                 succ[ pos + 1 ]:= result[ pos + 1 ] + 1;
                 VectorIterator_ResetPrefix( succ, iter!.m, pos,
                                         Sum( succ{ [ 1 .. pos ] } ) );
               fi;
               iter!.next:= succ;
               return result;
             end,

             IsDoneIterator:= iter -> ( iter!.next = false ),

             ShallowCopy:= function( iter )
             end,

             m:= m,
             n:= n,
             s:= s,

             next:= next,
           ) );
    end );

#############################################################################
##
#F  MinimalDegreeHardGAP( <q>, <n>, <e> )
##
##  - assume that the cheap criteria do not yield the result
##    (e.g., m = 2 cannot occur here)
##  - assume q < e
##  - assume that n equals OrderMod( q, e )
##  - do not use the cache
##  - do not create the monomials by computing all q-adic expansions at once,
##    but create these coefficient vectors one by one, using an iterator
##  - for increasing values of m, run over those coefficient vectors
##    of length n and with entries in \{ 0, 1, \ldots, m \}
##    such that the coefficient sum is s
##    and such that a largest entry is in the first position.
##
InstallGlobalFunction( MinimalDegreeHardGAP, function( q, n, e )
    local z, e1, m, powers, a, v;

    Assert( 2, PowerMod( q, n, e ) = 1 and OrderModExt( q, e, n ) = n );

    # If A[q,n,z] has small dimension and e is not very small
    # then enumerate its basis.
    # (For small e, we known that the m values are small,
    # so it is cheaper not to enumerate.)
    if e > 100 then
      z:= ( q^n - 1 ) / e;
      if z <= 10^5 then
        return LoewyStructureInfoGAP( q, n, z ).m;
      fi;
    fi;

    # m(q,e) is divisible by \gcd( e, q-1 ), and m(q,e) \geq 3.
    e1:= Gcd( e, q-1 );
    if e1 = 1 then
      m:= 3;
    elif e1 = 2 then
      m:= 4;
    else
      m:= e1;
    fi;

    powers:= List( [ 1 .. n-1 ], i -> PowerModInt( q, i, e ) );

    while true do
      for a in [ 1 .. m ] do
        for v in VectorIterator( a, n-1, m-a ) do
          if ( a + v * powers ) mod e = 0 then
            return m;
          fi;
        od;
      od;
      m:= m + e1;
    od;
    end );

#############################################################################
##
#M  MinimalDegreeOfSingerAlgebraJulia(  )
#M  MinimalDegreeOfSingerAlgebraJulia( <q>, <e> )
#M  MinimalDegreeOfSingerAlgebraJulia( <q>, <n>, <e> )
##
##  Use the same criteria as for the &GAP; variant.
##  These methods are available only if &Julia; is available.
##
if IsPackageMarkedForLoading( "JuliaInterface", "" ) then

InstallMethod( MinimalDegreeOfSingerAlgebraJulia,
    [ "IsSingerAlgebra and HasLoewyStructureInfoJulia" ],
    A -> JuliaToGAP( IsInt,
             Julia.Base.get( LoewyStructureInfoJulia( A ),
                             JuliaSymbol( "m" ), 0 )  ));

InstallMethod( MinimalDegreeOfSingerAlgebraJulia,
    [ "IsSingerAlgebra" ],
    function( A )
    local paras, q, n, e;

    # Note that we cannot assume that the parameter 'n' is equal to
    # the multiplicative order of 'q' modulo 'e'.
    paras:= ParametersOfSingerAlgebra( A );
    q:= paras[1];
    n:= paras[2];
    e:= SingerAlgE( A );

    # Supply 'A' as an argument in order to set the GAP attribute.
    return JuliaToGAP( IsInt,
               Julia.SingerAlg.MinimalDegree( q, OrderModExt( q, e, n ), e,
                                              A ) );
    end );

InstallMethod( MinimalDegreeOfSingerAlgebraJulia,
    [ "IsPosInt", "IsPosInt" ],
    { q, e } -> JuliaToGAP( IsInt,
                    Julia.SingerAlg.MinimalDegree( q, GAPToJulia( e ) ) ) );

InstallMethod( MinimalDegreeOfSingerAlgebraJulia,
    [ "IsPosInt", "IsPosInt", "IsPosInt" ],
    { q, n, e } -> JuliaToGAP( IsInt,
               Julia.SingerAlg.MinimalDegree( q, n, GAPToJulia( e ) ) ) );

fi;

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

max_stars_repo_name: oscar-system/SingerAlg

max_stars_repo_licenses: Naumen Condor-1.1 MS-PL

max_issues_repo_name: oscar-system/SingerAlg

max_issues_repo_licenses: Naumen Condor-1.1 MS-PL

max_forks_repo_name: oscar-system/SingerAlg

max_forks_repo_licenses: Naumen Condor-1.1 MS-PL


text:
n:=5^(4^(3^2));;
s := String(n);;
m := Length(s);
# 183231
s{[1..20]};
# "62060698786608744707"
s{[m-19..m]};
# "92256259918212890625"


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

@TInt := @.cond(x->x.t=TInt);
@TReal := @.cond(x->IsRealT(x.t));
@_scalar := @.cond(x->IsOrdT(x.t) or IsRealT(x.t) or ObjId(x.t)=TPtr);
@TVect := @.cond(x->IsVecT(x.t));
@Value := @.cond(x->IsValue(x));
@nth := @.cond(x->ObjId(x)=nth);

_isa := self -> self.opts.vector.isa;
_epi := (self, o) -> Concat("epi", self.ctype_suffixval(o.t, _isa(self)));
_px := (self, o) -> self.ctype_suffixval(o.t, _isa(self));
_vp := (o) -> let(pp := Last(o.args),
    Cond(ObjId(pp)=vparam, pp.p,
	 IsValue(pp), pp.v,
	 pp));

Class(AVXUnparser, SSEUnparser, rec(
    # --------------------------------
    # ISA constructs, general
    # -------------------------------

    # This is a general suffix for intrinsics that is determine from the data type
    ctype_suffix := (self, t, isa) >> Cond(
        t = TVect(T_Real(64), 4), "pd",
        t = TVect(T_Real(32), 8), "ps",
        t = TVect(TReal, 4) and isa=AVX_4x64f, "pd",
        t = TVect(TReal, 8) and isa=AVX_8x32f, "ps",
        Inherited(t, isa)
    ),

    ctype_prefix := (self, t) >> Cond( _avxT(t, self.opts), "_mm256", "_mm" ),

    # This is the type used for declarations of vector variables
    ctype := (self, t, isa) >> Cond(
        t in [TReal, TVect(TReal, 1)],
          Cond(
            isa = AVX_4x64f, "double",
            isa = AVX_8x32f, "float",
            "UNKNOWN_TYPE"),
        t = T_Real(64), "double",
        t = T_Real(32), "float",
        # else
          Cond(
            t = TVect(TReal, 2),
              Cond(
                isa = AVX_4x64f, "__m128d",
                isa = AVX_8x32f, "__m64",
                "UNKNOWN_TYPE"),
            t = TVect(TReal, 4),
              Cond(
                isa = AVX_4x64f, "__m256d",
                isa = AVX_8x32f, "__m128",
                "UNKNOWN_TYPE"),
            t = TVect(TReal, 8),
              Cond(
                isa = AVX_8x32f, "__m256",
                "UNKNOWN_TYPE"),
            t = TVect(T_Real(64), 4), "__m256d",
            t = TVect(T_Real(32), 8), "__m256",
            t = TVect(T_Int(32),  8), "__m256i",
            t = TVect(T_UInt(32), 8), "__m256i",
            t = TVect(T_Real(64), 2), "__m128d",
            t = TVect(T_Real(32), 4), "__m128",
            t = TVect(T_Real(32), 2), "__m64",
            Inherited(t, isa))
    ),

    cvalue_suffix  := (self, t)  >> let( isa := _isa(self), Cond(
        (t = TReal and isa = AVX_8x32f) or t = T_Real(32), "f",
        (t = TReal and isa = AVX_4x64f) or t = T_Real(64), "",
        Inherited(t)
    )),

    vhex := (self, o, i, is) >> Print("_mm_set_", _epi(self, o), "(", self.infix(Reversed(o.p), ", "), ")"),

    vparam := (self, o, i, is) >> When(Length(o.p)=1, Print(o.p[1]), iclshuffle(o.p)),

    Value := (self, o, i, is) >> Cond(
        o.t = TString, Print(o.v),

        o.t = TReal or ObjId(o.t) = T_Real, let(v := When(IsCyc(o.v), ReComplex(Complex(o.v)), Double(o.v)),
            When(v<0, Print("(", v, self.cvalue_suffix(o.t), ")"), Print(v, self.cvalue_suffix(o.t)))),

        o.t = TComplex, Print("COMPLEX(", ReComplex(Complex(o.v)), self.cvalue_suffix(TReal), ", ", ImComplex(Complex(o.v)), self.cvalue_suffix(TReal), ")"),

	o.t in [TInt, TUChar, TChar],
            When(o.v < 0, Print("(", o.v, ")"), Print(o.v)),

        ObjId(o.t) = TVect and Length(Set(o.v)) = 1,
            Cond( self.cx.isInside(Value) and Length(self.cx.Value) >= 2, # nested in an array
		    Print("{", self.infix(Replicate(o.t.size, o.v[1]), ", "), "}"),
                  # else
                    Cond( _avxT(o.t, self.opts), let( sfx := self.ctype_suffix(o.t, _isa(self)), pfx := self.ctype_prefix(o.t),
                        Cond( o.v[1] = 0,
                            self.printf("$1_setzero_$2()", [pfx, sfx]),
                            self.printf("$1_set1_$2($3)", [pfx, sfx, o.v[1]]))),
                        Inherited(o, i, is))),
        ObjId(o.t) = TVect,
            Cond( self.cx.isInside(Value) and Length(self.cx.Value) >= 2, # nested in an array
		    Print("{", self.infix((o.v), ", "), "}"),
                # else
                    Cond( _avxT(o.t, self.opts), let( sfx := self.ctype_suffix(o.t, _isa(self)), pfx := self.ctype_prefix(o.t),
		        Print(pfx, "_set_", sfx, "(", self.infix(Reversed(o.v), ", "), ")")),
                        Inherited(o, i, is))),

        IsArray(o.t),
           Print("{", self.infix(o.v, ", "), "}"),

        ObjId(o.t) = TSym,
            Print("(", self.declare(o.t, [], 0, 0), ") ", o.v),

	o.t = TBool, Print(When(o.v, "1", "0")),

	#Error(self,".Value cannot unparse type ",o.t)
        Inherited(o, i, is)
    ),

    vpack := (self, o, i, is) >> Cond( _avxT(o.t, self.opts),
            Print("_mm256_set_", self.ctype_suffix(o.t, _isa(self)), "(", self.infix(Reversed(o.args), ", "), ")"),

        Inherited(o, i, is)),

    vdup := (self, o, i, is) >> CondPat(o,
            [vdup, @(1, [nth, deref]), @TInt], let( isa := _isa(self),
                t := o.args[1].t, pfx := self.ctype_prefix(o.t),
                Cond( t = T_Real(32) or (t=TReal and isa=AVX_8x32f),
                        self.printf("$1_broadcast_ss($2)", [pfx, o.args[1].toPtr(t)]),
                      t = T_Real(64) or (t=TReal and isa=AVX_4x64f),
                        self.printf("$1_broadcast_sd($2)", [pfx, o.args[1].toPtr(t)]),
                      t = TVect(T_Real(32), 2) or (t=TVect(TReal, 2) and isa=AVX_8x32f),
                        self.printf("$1_castpd_ps($1_broadcast_sd($2))", [pfx, o.args[1].toPtr(T_Real(64))]),
                      t = TVect(T_Real(64), 2) or (t=TVect(TReal, 2) and isa=AVX_4x64f),
                        self.printf("$1_broadcast_pd($2)", [pfx, o.args[1].toPtr(t)]),
                      t = TVect(T_Real(32), 4) or (t=TVect(TReal, 4) and isa=AVX_8x32f),
                        self.printf("$1_castpd_ps($1_broadcast_pd($2))", [pfx, o.args[1].toPtr(TVect(T_Real(64), 2))]),
                      Error("unexpected vdup load"))),
            #[vdup, @nth,@.cond(x->x.t=TInt and x.v=2)],
            #    Print("_mm_loaddup_", sfx, "(&(", self(o.args[1], i, is), "))"),
            [vdup, @, @TInt],
                Print(When(_isa(self).v = 4, "_mm256_set1_pd", "_mm256_set1_ps"), "(", self(o.args[1], i, is), ")")
        ),

    # Declarations
    TVect := (self, t, vars, i, is) >> let( ctype := self.ctype(t, _isa(self)),
        Print(ctype, " ", self.infix(vars, ", "))),

    TReal := ~.TVect,

    TInt := (self, t, vars, i, is) >> Print("int ", self.infix(vars, ", ")),

    TBool := (self, t, vars, i, is) >> Print("BOOL ", self.infix(vars, ", ")),

    # Arithmetic
    mul := (self, o, i, is) >> let(n := Length(o.args), Cond(
        not IsVecT(o.t),
            Print("(", self.pinfix(o.args, ")*("), ")"),
        not _avxT(o.t, self.opts),
            Inherited(o, i, is),
	n > 2 and n mod 2 <> 0,
            self(mul(o.args[1], ApplyFunc(mul, Drop(o.args, 1))), i, is),
        n > 2,
            self(mul(ApplyFunc(mul, o.args{[1..n/2]}), ApplyFunc(mul, o.args{[n/2+1..n]})), i, is),
        let( sfx := self.ctype_suffix(o.t, _isa(self)),
         CondPat(o,
           [mul, @TReal, @TVect],
              self(mul(vdup(o.args[1],o.t.size), o.args[2]), i, is),
           [mul, @TVect, @TReal],
              self(mul(o.args[1], vdup(o.args[2],o.t.size)), i, is),
           [mul, @TInt, @TVect],
              self(mul(vdup(_toReal(o.args[1]),o.t.size), o.args[2]), i, is),
           [mul, @TVect, @TInt],
              self(mul(o.args[1], vdup(_toReal(o.args[2]),o.t.size)), i, is),
           [mul, @TVect,   @TVect],
              self.printf("_mm256_mul_$1($2, $3)", [sfx, o.args[1], o.args[2]]),
           Error("Don't know how to unparse <o>. Unrecognized type combination")))
    )),

    # -- add --
    add := (self, o, i, is) >> let(n := Length(o.args), Cond(
	not IsVecT(o.t),
            self.pinfix(o.args, " + "),
        not _avxT(o.t, self.opts),
            Inherited(o, i, is),
        n > 2 and n mod 2 <> 0,
            self(add(o.args[1], ApplyFunc(add, Drop(o.args, 1))), i, is),
        n > 2,
            self(add(ApplyFunc(add, o.args{[1..n/2]}), ApplyFunc(add, o.args{[n/2+1..n]})), i, is),
        let(sfx := self.ctype_suffix(o.t, _isa(self)), saturated:= When(_isa(self).isFixedPoint and _isa(self).saturatedArithmetic, "s", ""),
          CondPat(o,
            [add, @TReal, @TVect],
                self(add(vdup(o.args[1],o.t.size), o.args[2]), i, is),
            [add, @TVect, @TReal],
                self(add(o.args[1], vdup(o.args[2],o.t.size)), i, is),
            [add, @TInt, @TVect],
                self(add(vdup(_toReal(o.args[1]),o.t.size), o.args[2]), i, is),
            [add, @TVect, @TInt],
                self(add(o.args[1], vdup(_toReal(o.args[2]),o.t.size)), i, is),
            [add, @TVect,   @TVect],
                self.printf("_mm256_add$1_$2($3, $4)", [saturated, sfx, o.args[1], o.args[2]]),
            Error("Don't know how to unparse <o>. Unrecognized type combination")))
    )),

    sub := (self, o, i, is) >> Cond( _avxT(o.t, self.opts), let(
            isa       := _isa(self),
            sfx       := self.ctype_suffix(o.t, isa),
            saturated := When(isa.isFixedPoint and isa.saturatedArithmetic, "s", ""),
            self.printf("_mm256_sub$1_$2($3, $4)", [saturated, sfx, o.args[1], o.args[2]])),
        # else
            Inherited(o, i, is)),

    neg := (self, o, i, is) >> Cond( _avxT(o.t, self.opts),
        self(mul(neg(o.t.one()), o.args[1]), i, is),
        Inherited(o, i, is)),

    stickyNeg := ~.neg,

    sqrt  := (self, o, i, is) >> Cond( _avxT(o.t, self.opts),
        Checked( IsRealT(o.t.t), self.printf("_mm256_sqrt_$1($2)", [self.ctype_suffix(o.t, _isa(self)), o.args[1]])),
        Inherited(o, i, is)),

    rsqrt := (self, o, i, is) >> Cond( _avxT(o.t, self.opts),
        let( sfx := self.ctype_suffix(o.t, _isa(self)),
            Checked( sfx="ps", self.printf("_mm256_rsqrt_ps($1)", [o.args[1]]))),
        Inherited(o, i, is)),

    max := (self, o, i, is) >> let(n := Length(o.args), Cond(
        not _avxT(o.t, self.opts),
            Inherited(o, i, is),
	n > 2 and n mod 2 <> 0,
            self(max(o.args[1], ApplyFunc(max, Drop(o.args, 1))), i, is),
        n > 2,
            self(max(ApplyFunc(max, o.args{[1..n/2]}), ApplyFunc(max, o.args{[n/2+1..n]})), i, is),
        let( sfx := self.ctype_suffix(o.t, _isa(self)),
         CondPat(o,
           [max, @TReal, @TVect],
              self(max(vdup(o.args[1],o.t.size), o.args[2]), i, is),
           [max, @TVect, @TReal],
              self(max(o.args[1], vdup(o.args[2],o.t.size)), i, is),
           [max, @TInt, @TVect],
              self(max(vdup(_toReal(o.args[1]),o.t.size), o.args[2]), i, is),
           [max, @TVect, @TInt],
              self(max(o.args[1], vdup(_toReal(o.args[2]),o.t.size)), i, is),
           [max, @TVect,   @TVect],
              self.printf("_mm256_max_$1($2, $3)", [sfx, o.args[1], o.args[2]]),
           Error("Don't know how to unparse <o>. Unrecognized type combination")))
    )),

    min := (self, o, i, is) >> let(n := Length(o.args), Cond(
        not _avxT(o.t, self.opts),
            Inherited(o, i, is),
	n > 2 and n mod 2 <> 0,
            self(min(o.args[1], ApplyFunc(min, Drop(o.args, 1))), i, is),
        n > 2,
            self(min(ApplyFunc(min, o.args{[1..n/2]}), ApplyFunc(min, o.args{[n/2+1..n]})), i, is),
        let( sfx := self.ctype_suffix(o.t, _isa(self)),
         CondPat(o,
           [min, @TReal, @TVect],
              self(min(vdup(o.args[1],o.t.size), o.args[2]), i, is),
           [min, @TVect, @TReal],
              self(min(o.args[1], vdup(o.args[2],o.t.size)), i, is),
           [min, @TInt, @TVect],
              self(min(vdup(_toReal(o.args[1]),o.t.size), o.args[2]), i, is),
           [min, @TVect, @TInt],
              self(min(o.args[1], vdup(_toReal(o.args[2]),o.t.size)), i, is),
           [min, @TVect,   @TVect],
              self.printf("_mm256_min_$1($2, $3)", [sfx, o.args[1], o.args[2]]),
           Error("Don't know how to unparse <o>. Unrecognized type combination")))
    )),

    # assuming we have ICC <ia32intrin.h> here
    log := (self, o, i, is) >> Cond( _avxT(o.t, self.opts), let( sfx := self.ctype_suffix(o.t, _isa(self)),
        Checked( sfx in ["ps", "pd"], Cond(
            Length(o.args)>1 and (o.args[2]=2 or o.args[2]=o.t.value(2)),
                self.printf("_mm256_log2_$1($2)", [sfx, o.args[1]]),
            Length(o.args)>1 and (o.args[2]=10 or o.args[2]=o.t.value(10)),
                self.printf("_mm256_log10_$1($2)", [sfx, o.args[1]]),
            Length(o.args)=1 or o.args[2]=d_exp(1) or o.args[2]=o.t.value(d_exp(1)),
                self.printf("_mm256_log_$1($2)", [sfx, o.args[1]]),
            self.printf("_mm256_div_$1(_mm256_log_$1($2), _mm256_log_$1($3))", [sfx, o.args[1]])))),
        Inherited(o, i, is)),

    # assuming we have ICC <ia32intrin.h> here
    exp := (self, o, i, is) >> Cond( _avxT(o.t, self.opts), let( sfx := self.ctype_suffix(o.t, _isa(self)),
        Checked( sfx in ["ps", "pd"], self.printf("_mm256_exp_$1($2)", [sfx, o.args[1]]))),
        Inherited(o, i, is)),

    # assuming we have ICC <ia32intrin.h> here
    pow := (self, o, i, is) >> Cond( _avxT(o.t, self.opts), let( sfx := self.ctype_suffix(o.t, _isa(self)),
        Checked( sfx in ["ps", "pd"], Cond(
            o.args[1]=2 or o.args[1]=o.t.value(2),
                self.printf("_mm256_exp2_$1($2)", [sfx, o.args[2]]),
            o.args[1]=d_exp(1) or o.args[1]=o.t.value(d_exp(1)),
                self.printf("_mm256_exp_$1($2)", [sfx, o.args[2]]),
            self.printf("_mm256_pow_$1($2, $3)", [sfx, o.args[1], o.args[2]])))),
        Inherited(o, i, is)),

    # --------------------------------
    # logic
    #
    bin_xor := (self, o, i, is) >> Cond( _avxT(o.t, self.opts),
        self.printf("_mm256_xor_$1($2, $3)", [self.ctype_suffix(o.t, _isa(self)), o.args[1], o.args[2]]),
        Inherited(o, i, is)),
    bin_and := (self, o, i, is) >> Cond( _avxT(o.t, self.opts),
        self.printf("_mm256_and_$1($2, $3)", [self.ctype_suffix(o.t, _isa(self)), o.args[1], o.args[2]]),
        Inherited(o, i, is)),
    bin_or := (self, o, i, is) >> Cond( _avxT(o.t, self.opts),
        self.printf("_mm256_or_$1($2, $3)", [self.ctype_suffix(o.t, _isa(self)), o.args[1], o.args[2]]),
        Inherited(o, i, is)),
    bin_andnot := (self, o, i, is) >> Cond( _avxT(o.t, self.opts),
        self.printf("_mm256_andnot_$1($2, $3)", [self.ctype_suffix(o.t, _isa(self)), o.args[1], o.args[2]]),
        Inherited(o, i, is)),

    # --------------------------------
    # ISA specific : AVX_4x64f
    #
    cmpge_4x64f := (self, o, i, is) >> self.prefix("_mm256_cmp_pd", Concat(o.args, ["_CMP_GE_OQ"])),
    logic_and_4x64f := (self, o, i, is) >> self.prefix("_mm256_and_pd", o.args),
    logic_xor_4x64f := (self, o, i, is) >> self.prefix("_mm256_xor_pd", o.args),

    vloadu_4x64f   := (self, o, i, is) >> self.prefix("_mm256_loadu_pd", o.args),
    vstoreu_4x64f  := (self, o, i, is) >> Print(Blanks(i), self.prefix("_mm256_storeu_pd", o.args), ";\n"),

    vinsert_2l_4x64f := (self, o, i, is) >> self.prefix("_mm256_insertf128_pd", o.args),
    vloadmask_4x64f  := (self, o, i, is) >> self.printf("_mm256_maskload_pd($1, _mm256_set_epi32($2))", [o.args[1], ()->PrintDel(_vp(o), ", ")]),
    vbroadcast_4x64f := (self, o, i, is) >> self.printf("_mm256_broadcast_sd($1)", [o.args[1]]),
    # doublecheck this
    vextract_2l_4x64f := (self, o, i, is) >> self.prefix("_mm256_extractf128_pd", o.args),

    vstore_2l_4x64f := (self, o, i, is) >> self.prefix("_mm256_extractf128_pd", o.args),
    vstoremask_4x64f := (self, o, i, is) >> Print(Blanks(i), self.printf("_mm256_maskstore_pd($1, _mm256_set_epi32($3), $2)", [o.args[1], o.args[2], ()->PrintDel(List(_vp(o), e->e.v), ", ")]), ";\n"),

    vunpacklo_4x64f := (self, o, i, is) >> self.prefix("_mm256_unpacklo_pd", o.args),
    vunpackhi_4x64f := (self, o, i, is) >> self.prefix("_mm256_unpackhi_pd", o.args),
    vpermf128_4x64f := (self, o, i, is) >> self.printf("_mm256_permute2f128_pd($1, $2, ($3) | (($4) << 4))",
        let(l := _vp(o)-1, [o.args[1], o.args[2], l[1], l[2]])),
    vshuffle_4x64f := (self, o, i, is) >> self.printf("_mm256_shuffle_pd($1, $2, ($3) | (($4) << 1) | (($5) << 2) | (($6) << 3))",
        let(l := _vp(o)-1, [o.args[1], o.args[2], l[1], l[2], l[3], l[4]])),
    vperm_4x64f := (self, o, i, is) >> self.printf("_mm256_permute_pd($1, ($2) | (($3) << 1) | (($4) << 2) | (($5) << 3))",
        let(l := _vp(o)-1, [o.args[1], l[1], l[2], l[3], l[4]])),
    vblend_4x64f := (self, o, i, is) >> self.printf("_mm256_blend_pd($1, $2, ($3) | (($4) << 1) | (($5) << 2) | (($6) << 3))",
        let(l := _vp(o)-1, [o.args[1], o.args[2], l[1], l[2], l[3], l[4]])),

    vuunpacklo_4x64f := (self, o, i, is) >> self(o.toBinop(), i, is),
    vuunpackhi_4x64f := (self, o, i, is) >> self(o.toBinop(), i, is),
    vupermf128_4x64f := (self, o, i, is) >> self(o.toBinop(), i, is),
    vushuffle_4x64f  := (self, o, i, is) >> self(o.toBinop(), i, is),
    vuperm2_4x64f    := (self, o, i, is) >> self(o.toBinop(), i, is),

    addsub_4x64f := (self, o, i, is) >> When(Length(o.args) > 2,
        Error("addsub_2x64f is strictly binary"),
        CondPat(o,
           [addsub_4x64f, @TReal, @TVect],
              self(_computeExpType(addsub_4x64f(vdup(o.args[1],o.t.size), o.args[2])), i, is),
           [addsub_4x64f, @TVect, @TReal],
              self(_computeExpType(addsub_4x64f(o.args[1], vdup(o.args[2],o.t.size))), i, is),
           [addsub_4x64f, @TInt, @TVect],
              self(_computeExpType(addsub_4x64f(vdup(_toReal(o.args[1]),o.t.size), o.args[2])), i, is),
           [addsub_4x64f, @TVect, @TInt],
              self(_computeExpType(addsub_4x64f(o.args[1], vdup(_toReal(o.args[2]),o.t.size))), i, is),
           [addsub_4x64f, @TVect,   @TVect],
              self.printf("_mm256_addsub_pd($1, $2)", [o.args[1], o.args[2]]),
           Error("Don't know how to unparse <o>. Unrecognized type combination")
        )),

    fmaddsub_4x64f := (self, o, i, is) >> When(
        Length(o.args) <> 3, Error("fmaddsub_4x64f is strictly ternary"),
        CondPat(o,
           [fmaddsub_4x64f, @TVect, @TVect, @TVect],
              self.printf("_mm256_fmaddsub_pd($1, $2, $3, 0)", [o.args[1], o.args[2], o.args[3]]),
           Error("Don't know how to unparse <o>. Unrecognized type combination")
        )),

    vzero_4x64f := (self, o, i, is) >> Print("_mm256_setzero_pd()"),

    # --------------------------------
    # ISA specific : AVX_8x32f
    #
    logic_and_8x32f := (self, o, i, is) >> self.prefix("_mm256_and_ps", o.args),
    logic_xor_8x32f := (self, o, i, is) >> self.prefix("_mm256_xor_ps", o.args),

    vloadu_8x32f      := (self, o, i, is) >> self.prefix("_mm256_loadu_ps", o.args),
    vstoreu_8x32f     := (self, o, i, is) >> Print(Blanks(i), self.prefix("_mm256_storeu_ps", o.args), ";\n"),

    vinsert_4l_8x32f  := (self, o, i, is) >> self.prefix("_mm256_insertf128_ps", o.args),
    vloadmask_8x32f   := (self, o, i, is) >> self.printf("_mm256_maskload_ps($1, _mm256_set_epi32($2))", [o.args[1], ()->PrintDel(_vp(o), ", ")]),
    vhdup_8x32f       := (self, o, i, is) >> self.printf("_mm256_movehdup_ps($1)", [o.args[1]]),
    vldup_8x32f       := (self, o, i, is) >> self.printf("_mm256_moveldup_ps($1)", [o.args[1]]),

    vextract_4l_8x32f := (self, o, i, is) >> self.prefix("_mm256_extractf128_ps", o.args),
    vstore_4l_8x32f   := (self, o, i, is) >> self.prefix("_mm256_extractf128_ps", o.args),
    vstoremask_8x32f  := (self, o, i, is) >> Print(Blanks(i), self.printf("_mm256_maskstore_ps($1, _mm256_set_epi32($3), $2)", [o.args[1], o.args[2], ()->PrintDel(List(_vp(o), e->e.v), ", ")]), ";\n"),

    vunpacklo_8x32f   := (self, o, i, is) >> self.prefix("_mm256_unpacklo_ps", o.args),
    vunpackhi_8x32f   := (self, o, i, is) >> self.prefix("_mm256_unpackhi_ps", o.args),

    vpermf128_8x32f   := (self, o, i, is) >> self.printf("_mm256_permute2f128_ps($1, $2, ($3) | (($4) << 4))",
        let(l := _vp(o)-1, [o.args[1], o.args[2], l[1], l[2]])),

    vshuffle_8x32f    := (self, o, i, is) >> self.printf("_mm256_shuffle_ps($1, $2, ($3) | (($4) << 2) | (($5) << 4) | (($6) << 6))",
        let(l := _vp(o)-1, [o.args[1], o.args[2], l[1], l[2], l[3], l[4]])),

    vperm_8x32f       := (self, o, i, is) >> self.printf("_mm256_permute_ps($1, ($2) | (($3) << 2) | (($4) << 4) | (($5) << 6))",
        let(l := _vp(o)-1, [o.args[1], l[1], l[2], l[3], l[4]])),
    vblend_8x32f := (self, o, i, is) >> self.printf("_mm256_blend_ps($1, $2, ($3) | (($4) << 1) | (($5) << 2) | (($6) << 3) | (($7) << 4) | (($8) << 5) | (($9) << 6) | (($10) << 7))",
        let(l := _vp(o)-1, [o.args[1], o.args[2], l[1], l[2], l[3], l[4], l[5], l[6], l[7], l[8]])),

    vuunpacklo_8x32f := (self, o, i, is) >> self(o.toBinop(), i, is),
    vuunpackhi_8x32f := (self, o, i, is) >> self(o.toBinop(), i, is),
    vupermf128_8x32f := (self, o, i, is) >> self(o.toBinop(), i, is),
    vushuffle_8x32f  := (self, o, i, is) >> self(o.toBinop(), i, is),

    addsub_8x32f := (self, o, i, is) >> When(Length(o.args) > 2,
        Error("addsub_2x64f is strictly binary"),
        CondPat(o,
           [addsub_8x32f, @TReal, @TVect],
              self(_computeExpType(addsub_8x32f(vdup(o.args[1],o.t.size), o.args[2])), i, is),
           [addsub_8x32f, @TVect, @TReal],
              self(_computeExpType(addsub_8x32f(o.args[1], vdup(o.args[2],o.t.size))), i, is),
           [addsub_8x32f, @TInt, @TVect],
              self(_computeExpType(addsub_8x32f(vdup(_toReal(o.args[1]),o.t.size), o.args[2])), i, is),
           [addsub_8x32f, @TVect, @TInt],
              self(_computeExpType(addsub_8x32f(o.args[1], vdup(_toReal(o.args[2]),o.t.size))), i, is),
           [addsub_8x32f, @TVect,   @TVect],
              self.printf("_mm256_addsub_ps($1, $2)", [o.args[1], o.args[2]]),
           Error("Don't know how to unparse <o>. Unrecognized type combination")
        )),

    fmaddsub_8x32f := (self, o, i, is) >> When(
        Length(o.args) <> 3, Error("fmaddsub_8x32f is strictly ternary"),
        CondPat(o,
           [fmaddsub_8x32f, @TVect, @TVect, @TVect],
              self.printf("_mm256_fmaddsub_ps($1, $2, $3, 0)", [o.args[1], o.args[2], o.args[3]]),
           Error("Don't know how to unparse <o>. Unrecognized type combination")
        )),

    vzero_8x32f := (self, o, i, is) >> Print("_mm256_setzero_ps()"),

    #--------------------------------
    # Conversion

    vcvt_8x32f_4x64f  := (self, o, i, is) >> self.printf("_mm256_cvtpd_ps($1)", [o.args[1]]),
    vcvt_4x64f_4x32f  := (self, o, i, is) >> self.printf("_mm256_cvtps_pd($1)", [o.args[1]]),
    vcvt_4x64f_4x32i  := (self, o, i, is) >> self.printf("_mm256_cvtepi32_pd($1)", [o.args[1]]),
    vcvt_4x32i_4x64f  := (self, o, i, is) >> self.printf("_mm256_cvtpd_epi32($1)", [o.args[1]]),
    vcvtt_4x32i_4x64f := (self, o, i, is) >> self.printf("_mm256_cvttpd_epi32($1)", [o.args[1]]),
    vcvt_8x32f_8x32i  := (self, o, i, is) >> self.printf("_mm256_cvtepi32_ps($1)", [o.args[1]]),
    vcvt_8x32i_8x32f  := (self, o, i, is) >> self.printf("_mm256_cvtps_epi32($1)", [o.args[1]]),
    vcvtt_8x32i_8x32f := (self, o, i, is) >> self.printf("_mm256_cvttps_epi32($1)", [o.args[1]]),

));


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

Class(RulesCxRC_Op, RuleSet);
RewriteRules(RulesCxRC_Op, rec(

    vRC_Compose := Rule([vRC, @(1, Compose)], e -> Compose(List(@(1).val.children(), vRC))),
    vRC_SUM := Rule([vRC, @(1, SUM)], e -> SUM(List(@(1).val.children(), vRC))),
    vRC_SUMAcc := Rule([vRC, @(1, SUMAcc)], e -> SUMAcc(List(@(1).val.children(), vRC))),

    vRC_Container := Rule([vRC, @(1, [BB,Buf,Inplace,Grp,NoPull,NoPullLeft,NoPullRight, NoDiagPullin, NoDiagPullinLeft, NoDiagPullinRight ])],
        e -> ObjId(@(1).val)(vRC(@(1).val.child(1)))),

    vRC_Data := Rule([vRC, @(1, Data)], e -> Data(@(1).val.var, @(1).val.value, vRC(@(1).val.child(1)))),

    vRC_RStep := Rule([vRC, @(1, RecursStep)], e -> RecursStep(2*@(1).val.yofs, 2*@(1).val.xofs,
                                                      vRC(@(1).val.child(1)))),

    vRC_ISum := Rule([vRC, @(1, ISum)], e -> ISum(@(1).val.var, @(1).val.domain, vRC(@(1).val.child(1)))),
    vRC_ICompose := Rule([vRC, @(1, ICompose)], e -> ICompose(@(1).val.var, @(1).val.domain, vRC(@(1).val.child(1)))),

    vRC_Grp := Rule([vRC, @(1, Grp)], e -> Grp(vRC(@(1).val.child(1)))),

    vRC_Scale := Rule([vRC, @(1, Scale)], e ->
         vRC(Diag(fConst(Rows(@(1).val), @(1).val.scalar))) * vRC(@(1).val.child(1))),

    vRC_CR := Rule([vRC, @(1, CR)], e -> @(1).val.child(1)),

    VTensor_VScale := Rule([@(1, VTensor), [VScale, @(2), @(3), @(4)]], e -> VScale(VTensor(@(2).val, @(1).val.vlen), @(3).val, @(4).val*@(1).val.vlen)),
    VTensor_VTensor := Rule([@(1, VTensor), @(2, VTensor)], e -> VTensor(@(2).val.child(1), @(1).val.vlen*@(2).val.vlen)),

    vRC_SMP := Rule([@(1, vRC), @(2, [SMPSum, SMPBarrier])],
        e -> let(s := @(2).val, CopyFields(s, rec(_children := List(s.children(), c->ObjId(@(1).val)(c)), dimensions := @(1).val.dimensions)))),

    # should RC and VRC be in some other place?
    vRC_TCvt := Rule( [@(0, [vRC, RC, VRC]), @(1, TCvt)],
        e -> let( t := @(1).val, TCvt( 2*t.n(), t.isa_to(), t.isa_from(), t.props()).withTags(t.getTags()).takeAobj(t) )),
));

Class(RulesCxRC_Term, RuleSet);
RewriteRules(RulesCxRC_Term, rec(
    vRC_VContainer := Rule([vRC, @(1, VContainer)], e -> VContainer(vRC(@(1).val.child(1)), @(1).val.isa.uncplx())),

#    RC_VTensor := Rule([ RC, @(1, VTensor) ], e -> RCVTensor(@(1).val.child(1), @(1).val.vlen)),
#
#    RCVTensor_VTensor := Rule(@(1, RCVTensor, e -> e.isReal()), e -> VTensor(@(1).val.child(1), 2*@(1).val.vlen)),

#
#    RC_VGathScat_sv := Rule([@(1, RC), @(2, [VGath_sv, VScat_sv])], e->
#        @(2).val.rcVariant(@(2).val.func, 2*@(2).val.v, @(2).val.sv)),
#
#    RC_VGathScat := Rule([@(1, RC), @(2, [VGath, VScat])], e->
#        ObjId(@(2).val)(@(2).val.func, 2*@(2).val.v)),
#
#    RC_VPrm_x_I := Rule([@(1, RC), @(2, VPrm_x_I)], e ->  VPrm_x_I(@(2).val.func, @(2).val.inverse, 2*@(2).val.v)),
#
#    RC_VDiag_x_I := Rule([@(1, RC), @(2, VDiag_x_I)], e ->  VDiag_x_I(@(2).val.element, 2*@(2).val.v)),
#
#    VContainer_RCDiag := Rule(@@(1, RCDiag, (e, cx) -> IsBound(cx.VContainer) and Length(cx.VContainer)=1),
#        (e, cx) -> RCVDiag(@@(1).val.element, Last(cx.VContainer).isa.getV())),
#
#    RCVTensor_Blk := Rule([@(1, RCVtensor), @(2, Blk)], e -> RCVDiag(RCData(@(2).val.element), 2*@(2).val.v)),

##--- RC rules moved to vRC

## Shared Franz/Yevgen Diag rules
    vRC_Diag      := Rule([@(1, vRC), @(2, Diag) ], e -> vRC(VDiag(@(2).val.element, 1))),

#############################
## Original Franz's rules:

    vRC_VDiag     := Rule([@(1, vRC), @(2, VDiag)], e -> let(v:=2*@(2).val.v,
        RCVDiag(VData(RCData(@(2).val.element), v), v))),

    vRC_VDiag_x_I := Rule([@(1, vRC), @(2, VDiag_x_I)], e -> let(
        v   := @(2).val.v,
        elt := @(2).val.element,
        dup := Cond(IsVecT(elt.range()), elt, VDup(elt, v)),
        RCVDiag(RCVData(dup), 2*v)
    )),

############################################
## Yevgen's searchable rules for autolib
    vRC_VDiag     := Rule([@(1, vRC), @(2, VDiag)], (e,cx) -> let(
	v   := 2*@(2).val.v,
	alt := Global.autolib.RCOND("diag", Global.autolib.rcondModeFreedom, maybe(),
		     RCVDiag(VData(RCData(@(2).val.element), v), v),
		     RCVDiagSplit(RCVDataSplit(@(2).val.element, v), v)),
	Cond( not IsBound(cx.opts.hack_vRC_VDiag), alt.child(2),
	      cx.opts.hack_vRC_VDiag = "search",   alt,
	      cx.opts.hack_vRC_VDiag = "compact",  alt.child(1),
	      cx.opts.hack_vRC_VDiag = "split",    alt.child(2) ))), # strip alternator outside of autolib

    vRC_VDiag_x_I := Rule([@(1, vRC), @(2, VDiag_x_I)], (e,cx) -> let(
        v   := @(2).val.v,
        elt := @(2).val.element,
	t   := elt.range(),
        # Cond(IsVecT(t), VDup(elt, 2), VDup(elt, 2*v)), v),  # does not work, dimension error
	Cond(IsRealT(t),
	         VDiag_x_I(elt, 2*v),

	     IsVecT(t) and IsRealT(t.t),
	         VDiag_x_I(VDup(elt, 2), 2*v), # NOTE: make sure sub-vector vdup is properly
		                               #        handled by the compiler rewrite rules
	     # else, complex
	     let(
		 dup := Cond(IsVecT(t), elt, VDup(elt, v)),
		 alt := Global.autolib.RCOND("diag", Global.autolib.rcondModeFreedom, maybe(),
                     RCVDiag(RCVData(dup), 2*v),
		     RCVDiagSplit(RCVDataSplitVec(dup), 2*v)),
		 Cond( not IsBound(cx.opts.hack_vRC_VDiag_x_I),  alt.child(2),  # strip alternator outside of autolib
		     cx.opts.hack_vRC_VDiag_x_I = "search",    alt,
		     cx.opts.hack_vRC_VDiag_x_I = "compact",   alt.child(1),
		     cx.opts.hack_vRC_VDiag_x_I = "split",     alt.child(2) )
	      )
	))),

    # we need to check the divisibility of dimensions, because all VGathXXX constructs output full vectors
    # This means that gather of 3 points won't be handled here.
    vRC_Gath := Rule([vRC, @@(1, Gath, (e,cx)->IsBound(cx.VContainer) and cx.VContainer<>[] and imod(2*e.func.domain(), Last(cx.VContainer).isa.getV())=0)],
        (e,cx) -> VGath_sv(@@(1).val.func, Last(cx.VContainer).isa.getV(), 2)),

    vRC_Scat := Rule([vRC, @@(1, [Scat,ScatAcc], (e,cx)->IsBound(cx.VContainer) and cx.VContainer<>[] and imod(2*e.func.domain(), Last(cx.VContainer).isa.getV())=0)],
        (e,cx) -> Cond(@@(1).val _is Scat,
	    VScat_sv(@@(1).val.func, Last(cx.VContainer).isa.getV(), 2),
	    VScat_svAcc(@@(1).val.func, Last(cx.VContainer).isa.getV(), 2))),

#    NOTE: why is this commented out?
#    VA: btw - in fixed size code this happens anyway with help of vRC_REALMAT and RulesVec through VTensor
#    vRC_Prm := Rule([vRC, @(1, Prm)], e -> VPrm_x_I(@(1).val.func, 2)),

    # NOTE: This rule looks evil, as explained below.
    #        Rules for VTensor( VRCDiag | VPRM | ... ) do not exist
    #        and if it fires prematurely (maybe due to some composite operator)
    #        the rewrite system will get stuck later.
    #
    #        This is exactly what happens in autolib ASP.DFT(N) using _CT1Tr.
    #
    #        Thus I blocked it. I think it should be possible to get rid of
    #        it in a graceful way.
    #
    vRC_REALMAT := Rule([vRC, @@(1).cond(
        (e,cx) -> not IsBound(cx.opts.autolib) and e.isReal() and not ObjId(e) in
            [VRCDiag, VPrm_x_I, VIxJ2, VJxI, VBlk,
             VGath_zero, VScat_zero, VGath_pc, VScat_pc, VScat_pcAcc, IxVGath_pc, IxVScat_pc,
             VGath_sv, VScat_sv, VScat_svAcc, VGath, VScat, VScatAcc, VGath_u, VScat_u, VScatAcc_u, VStretchGath, VStretchScat,
             Gath, Scat, ScatAcc])],
        e -> VTensor(@@(1).val, 2)),

    vRC_Blk := Rule([vRC, @(1, Blk, e -> not e.isReal())], e -> RCVBlk(MapMat(@(1).val.element, i->[i]), 2)),
    vRC_Blk1 := Rule([vRC, @(1, Blk1, e -> not e.isReal())], e -> vRC(Blk([[@(1).val.element]]))),

    VTensor_XXX := Rule([@(1,VTensor,e->e.vlen=2), @(2,[VGath, VScat, VScatAcc, VGath_sv, VScat_sv, VScat_svAcc, VGath_u, VScat_u, VScatAcc_u, VPrm_x_I])], e->vRC(e.child(1))),

    VTensor_I := Rule([@(1, VTensor), [I, @(2)]], e -> VPrm_x_I(fId(@(2).val), e.vlen)),

    vRC_VGath := Rule([@(1, vRC), @(2, VGath)], e -> VGath(@(2).val.func, 2*@(2).val.v)),
    vRC_VScat := Rule([@(1, vRC), @(2, [VScat, VScatAcc])], e -> ObjId(@(2).val)(@(2).val.func, 2*@(2).val.v)),

    vRC_VGath_u := Rule([@(1, vRC), @(2, VGath_u)], e -> VGath_u(fTensor(@(2).val.func, fBase(2, 0)), 2*@(2).val.v)),
    vRC_VScat_u := Rule([@(1, vRC), @(2, [VScat_u, VScatAcc_u])], e -> ObjId(@(2).val)(fTensor(@(2).val.func, fBase(2, 0)), 2*@(2).val.v)),

    vRC_VGath_sv := Rule([@(1, vRC), @(2, VGath_sv)], e -> VGath_sv(@(2).val.func, 2*@(2).val.v, 2*@(2).val.sv, Cond(@(2).val.rem _is Unk, Unk(TInt), 2*@(2).val.rem))),
    vRC_VScat_sv := Rule([@(1, vRC), @(2, [VScat_sv, VScat_svAcc])], e -> ObjId(@(2).val)(@(2).val.func, 2*@(2).val.v, 2*@(2).val.sv, Cond(@(2).val.rem _is Unk, Unk(TInt), 2*@(2).val.rem))),

    vRC_VStretchGath := Rule([@(1, vRC), @(2, VStretchGath)], e -> vRCStretchGath(@(2).val.func, @(2).val.part, 2*@(2).val.v)),
    vRC_VStretchScat := Rule([@(1, vRC), @(2, VStretchScat)], e -> vRCStretchScat(@(2).val.func, @(2).val.part, 2*@(2).val.v)),

    fix_VGathScat := Rule([@(1, [VGath, VScat, VScatAcc]), [@(2,fTensor), ...,
            [fId, @@(3).cond((e, cx)->IsBound(cx.VContainer) and cx.VContainer <> [] and @(1).val.v < Last(cx.VContainer).isa.getV() and IsInt(@(1).val.v*e / Last(cx.VContainer).isa.getV())
            and (not IsBound(cx.VTensor) or (IsBound(cx.VTensor) and cx.VTensor = []))
            and (not IsBound(cx.vRC) or (IsBound(cx.vRC) and cx.vRC = []))
	    and not (IsBound(cx.opts.assumeOuter_vRC) and cx.opts.assumeOuter_vRC)
             ) ]]],
        (e, cx) -> ObjId(@(1).val)(fTensor(DropLast(@(2).val.children(), 1), fId(@(1).val.v*@@(3).val / Last(cx.VContainer).isa.getV())), Last(cx.VContainer).isa.getV())),

    diagTensor_fPrecompute := Rule([diagTensor, [fPrecompute, @(1)], @(2)], x->fPrecompute(diagTensor(@(1).val, @(2).val))),
    # shouldn't this be in some other place:
    fStretch_fPrecompute := Rule([fStretch, @(1, fPrecompute), @(2), @(3)], e -> fPrecompute(fStretch(@(1).val.child(1), @(2).val, @(3).val))),

    RCVData_fPrecompute := Rule([@(0,[RCVData,RCVDataSplitVec,VDataRDup]), [fPrecompute, @(1)]],
	x->fPrecompute(ObjId(@(0).val)(@(1).val))),
    RCVData2_fPrecompute := Rule([@(0,[RCVDataSplit]), [fPrecompute, @(1)], @(2)],
	x->fPrecompute(ObjId(@(0).val)(@(1).val, @(2).val))),

    RCVData_fCompose    := Rule([RCVData, @(1, fCompose, x -> not IsType(Last(x.children()).range()))],
        e -> let(ch := @(1).val.children(),
	         chsplit := SplitBy(ch, x->IsType(x.range())),
		 fCompose(RCVData(fCompose(chsplit[1])), fCompose(chsplit[2])))),

    RCVDataSplitVec_fCompose    := Rule([RCVDataSplitVec, @(1, fCompose, x -> not IsType(Last(x.children()).range()))],
        e -> let(ch := @(1).val.children(),
	         chsplit := SplitBy(ch, x->IsType(x.range())),
		 fCompose(RCVDataSplitVec(fCompose(chsplit[1])), fTensor(fCompose(chsplit[2]), fId(2))))),

    # NOTE: potential problem -- RCVDataSplit_fCompose rule is missing

    RCVData_diagDirsum  := Rule([@(0,[RCVData,RCVDataSplitVec]), @(1,diagDirsum)],
        e -> ApplyFunc(diagDirsum, List(@(1).val.children(), ObjId(@(0).val)))),

    RCVDataSplit_diagDirsum  := Rule([ @(0, RCVDataSplit, x -> x.v mod 2 = 0),
                                   @(1, diagDirsum, x -> x.subdomainsDivisibleBy(@(0).val.v/2)),
                                   @(2)],
        e -> ApplyFunc(diagDirsum, List(@(1).val.children(), x->ObjId(@(0).val)(x, @(2).val)))),
    # pull out fPrecompute because it's stuck in there, RCVDataSplit cannot go through diagDirsum
    RCVDataSplit_diagDirsum_fPrecompute  := Rule( [@(1, diagDirsum), ..., @@(2, fPrecompute, (x, cx) -> cx.isInside(RCVDataSplit) and not @(1).val.subdomainsDivisibleBy(Last(cx.RCVDataSplit).v/2) ), ...],
        e -> fPrecompute(SubstTopDown(e, fPrecompute, ee -> ee.child(1)))),

    # When vector length is 2 we can go through compose
    RCVData2_fCompose_diagDirsum := Rule([RCVDataSplit, [@(1, fCompose), @(2, diagDirsum), ...], @(3).cond( x -> x = 2)],
        e -> fCompose(RCVDataSplit(@(2).val, 2), fTensor(fCompose(Drop(@(1).val.children(),1)), fId(2)))),

    VDup_fCompose := Rule([VDup, [@(1,fCompose), @(2), ...], @(3)],
        e -> fCompose(VDup(@(2).val, @(3).val), fCompose(Drop(@(1).val.children(),1)))),

    VDup_diagDirsum := Rule([VDup, @(1,diagDirsum), @(2)],
        e -> ApplyFunc(diagDirsum, List(@(1).val.children(), x->VDup(x, @(2).val)))),

    # YSV: Previously RCVBlk(...) was created with complex vector length,
    #      now we use real vlen for consistency
    vRC_VTensor_Blk := Rule([@(0, vRC), [@(1, VTensor), @(2, Blk)]],
        e -> let(v := @(1).val.vlen,
	         RCVBlk(MapMat(@(2).val.element, e -> Replicate(v, e)), 2*v))),

    vRC_VGath_pc := Rule([@(1, vRC), @(2, VGath_pc)],
        e -> let(g := @(2).val, VGath_pc(g.N*2, g.n*2, g.ofs*2, g.v*2))),

    vRC_VScat_pc := Rule([@(1, vRC), @(2, [VScat_pc, VScat_pcAcc])],
        e -> let(s := @(2).val, ObjId(@(2).val)(s.N*2, s.n*2, s.ofs*2, s.v*2))),

    vRC_IxVGath_pc := Rule([@(1, vRC), @(2, IxVGath_pc)],
        e -> let(g := @(2).val, IxVGath_pc(g.k, g.N*2, g.n*2, g.ofs*2, g.v*2))),

    vRC_IxVScat_pc := Rule([@(1, vRC), @(2, IxVScat_pc)],
        e -> let(s := @(2).val, IxVScat_pc(s.k, s.N*2, s.n*2, s.ofs*2, s.v*2))),

    vRC_VIxJ2 := Rule([vRC, @(1, VIxJ2)], e -> RCVIxJ2(@(1).val.v*2)),
    vRC_VJxI  := Rule([vRC, @(1, VJxI)], e -> VJxI(@(1).val.m, @(1).val.v*2)),

    vRC_VBlk_real := Rule([vRC, @(1, VBlk, e -> e.isReal())],
	e -> VBlk(MapMat(@(1).val.element, i ->
		TVect(i.t.t.realType(), 2*@(1).val.v).value(
		    ConcatList(i.v, j -> [j, j]))), 2*@(1).val.v)),

    vRC_VBlk_complex := Rule([vRC, @(1, VBlk, e -> not e.isReal())],
	e -> RCVBlk(MapMat(@(1).val.element, _unwrap), 2*@(1).val.v)),

    vRC_VGath_zero := Rule([vRC, @(1, VGath_zero)], e -> VGath_zero(@(1).val.N, @(1).val.n, @(1).val.v*2)),

    vRC_VScat_zero := Rule([vRC, @(1, VScat_zero)], e -> VScat_zero(@(1).val.N, @(1).val.n, @(1).val.v*2)),

    vRC_VPrm_x_I := Rule([vRC, @(1, VPrm_x_I)], e -> VPrm_x_I(@(1).val.func, @(1).val.v*2)),

    sv2__Scat_sv_fTensor := Rule([@(1, [VScat_sv, VScat_svAcc]), [@(2,fTensor), ..., [fId, @(3).cond(e -> IsInt(@(1).val.sv*EvalScalar(e)/@(1).val.v))]]],
        e -> let(v := @(1).val.v, sv := @(1).val.sv, n := EvalScalar(@(3).val), newn := sv*n/v,
	         scat := Cond(@(1).val _is VScat_sv, VScat, VScatAcc),
             scat(fTensor(DropLast(@(2).val.children(), 1), fId(newn)), v))),

    sv2__Gath_sv_fTensor := Rule([@(1, VGath_sv), [@(2,fTensor), ..., [fId, @(3).cond(e -> IsInt(@(1).val.sv*EvalScalar(e)/@(1).val.v))]]],
        e -> let(v := @(1).val.v, sv := @(1).val.sv, n := EvalScalar(@(3).val), newn := sv*n/v,
             VGath(fTensor(DropLast(@(2).val.children(), 1), fId(newn)), v))),

#    vRC_VRCDiag := Rule([vRC, [@(1, VRCDiag), @(2, VDup)]],
#        e -> let(v := @(1).val.v*@(2).val.v, VRCDiag(VDup(@(2).val.func, v), v))),
#
#    vRC_VRCDiag_fPrecompute := Rule([vRC, [@(1, VRCDiag), [fPrecompute, @(2, VDup)]]],
#        e -> let(v := @(1).val.v*@(2).val.v, VRCDiag(fPrecompute(VDup(@(2).val.func, v)), v)))
#

    # NOTE: this is a recipe for disaster
#    vRC_VRCDiag_VDup := Rule([vRC, @(1, VRCDiag, e->Length(Collect(e.element, VDup))=1)],
#        e -> let(vdup := Collect(@(1).val.element, VDup)[1], v := 2*vdup.v, VRCDiag(SubstTopDown(@(1).val.element, @(2, VDup), i->VDup(@(2).val.func, v)), v))),

    vRC_VRCDiag := Rule([vRC, @(1, VRCDiag)], e -> VRCDiag(VDataRDup(@(1).val.element), 2*@(1).val.v)),

    vRC_I := Rule( [@(0, vRC), @(1, I)],
        e -> I(2*@(1).val.params[1])),

    vRC_VOLMultiplication := Rule([vRC, @(1, VOLMultiplication)],
        e -> let(p := @(1).val.rChildren(), RCVOLMultiplication(p[1], p[2], 2*p[3]))),
    vRC_VOLConjMultiplication := Rule([vRC, @(1, VOLConjMultiplication)],
        e -> let(p := @(1).val.rChildren(), RCVOLConjMultiplication(p[1], p[2], 2*p[3]))),

));

RulesCxRC := MergedRuleSet(RulesCxRC_Op, RulesCxRC_Term);


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

_prepVars := function(vars, vargen, pfx)
    local v, res, newpfx, newv;
    res := []; newpfx := Concat("p", pfx);
    for v in vars do
        if IsArrayT(v.t) then
	    newv := vargen.next(TPtr(v.t.t), newpfx);
	    Add(res, [ v, newv, zallocate(newv, v.t) ]);
	else
	    Add(res, [ v, v, skip() ]);
	fi;
    od;
    return res;
end;

GenerateMultiIOBench := function(t,code, domain, range, opts)
    local vargen, code, inp, out, inout, args, allocs;
    vargen := VarGenNumeric();
    inp := List(domain, t -> vargen.next(t, "_in"));
    inp := _prepVars(inp, vargen, "_in");

    out := List(range, t -> vargen.next(t, "_out"));
    out := _prepVars(out, vargen, "_out");

    inout := Concatenation(out, inp);
    args := List(inout, x->x[2]);
    allocs := List(inout, x->x[3]);

    code := program(
	decl(args,
	    chain(
		SubstTopDown(code, program, x->chain(x.cmds)),
		func(TVoid, "init", [],
		    chain(allocs,
		    ApplyFunc(call, [var(opts.subInitName)]))),
		func(TVoid, "transform", [When(IsBound(opts.subParams), opts.subParams, [])],
		    ApplyFunc(call, Concatenation([var(opts.subName)], When(IsBound(opts.subParams), opts.subParams, []),args))))));
    return code;
end;

# Generate Bench does allocation for a program (a function and its initialization).
# It indirects both the init and the transform, does memory allocation for parameters
# in the new init and uses them to call the function in the new body.
#
#     init(){...}
#     transform(...){...}
#
# will be thus be renamed to:
#
#      <opts.subInitName>(){...}
#      <opts.subName>(...){...}
#
# and the following functions will be added:
#
#      init(){params=malloc();<opts.subInitName>();}
#      transform(){<opts.subName>(params);}
#
# GenerateBench(<t>,<code>,<opts>)
# <t> is the sums corresponding to the program <code> compiled with options <opts>
GenerateBench := (t,code,opts) ->
    GenerateMultiIOBench(t, code, t.dmn(), t.rng(), opts);


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

ObjLookup := arg -> ApplyFunc(ObjHash.objLookup, arg);
ObjAdd    := arg -> ApplyFunc(ObjHash.objAdd, arg);
MemClass  := arg -> ApplyFunc(ObjHash.memClass, arg);
MemClassFunc  := arg -> ApplyFunc(ObjHash.memClassFunc, arg);
SingletonAdd  := arg -> ApplyFunc(ObjHash.singletonAdd, arg);

objs := x -> Filtered(Collect(x, @), IsRec);
uids := x -> List(objs(x), o->o.uid);
### init
SingletonAdd(TInt);
SingletonAdd(TComplex);
SingletonAdd(TUnknown);
SingletonAdd(TReal);

MemClass(TArray);
MemClass(TArrayBase);
MemClass(TVect);

#ClassSPL.hash := ObjHash;
#Function.hash := ObjHash;

########
MemClass(Exp);
MemClass(ListableExp);
MemClass(Lambda);
MemClass(FList);
MemClass(FData);
MemClass(nth);
MemClass(Value);

MemClassFunc(Value, "new", "new_no_memo");
MemClassFunc(Value, "newbase", "newbase_no_memo");

max_stars_repo_name: franzfranchetti/spiral-software

max_issues_repo_name: franzfranchetti/spiral-software

max_forks_repo_name: franzfranchetti/spiral-software