text:
#############################################################################
##
#W  smlinfo.gi               GAP group library             Hans Ulrich Besche
##                                               Bettina Eick, Eamonn O'Brien
##
##  This file contains the ...
##

#############################################################################
##
#F  SMALL_GROUPS_INFORMATION
##
##  ...
SMALL_GROUPS_INFORMATION := [ ];

#############################################################################
##
#F  SmallGroupsInformation( size )
##
##  ...
InstallGlobalFunction( SmallGroupsInformation, function( size )
    local smav, idav, num, lib, t;

    if not IsPosInt(size) then
      ErrorNoReturn("usage: SmallGroupsInformation( size )");
    fi;

    smav := SMALL_AVAILABLE( size );
    idav := ID_AVAILABLE( size );

    if size = 1024 then
        Print( "The groups of size 1024 are not available. \n");
        return;
    fi;

    if smav = fail then
        Print( "The groups of size ", size, " are not available. \n");
        return;
    fi;

    lib := 1;
    if IsBound( smav.lib ) then
        lib := smav.lib;
    fi;

    if IsBound( smav.number ) then
        num := smav.number;
    else
        num := NUMBER_SMALL_GROUPS_FUNCS[ smav.func ]( size, smav ).number;
    fi;
    if num = 1 then
        Print("\n  There is 1 group of order ",size,".\n");
    else
        Print("\n  There are ",num," groups of order ",size,".\n" );
    fi;

    SMALL_GROUPS_INFORMATION[ smav.func ]( size, smav, num );

    Print("\n  This size belongs to layer ",lib,
          " of the SmallGroups library. \n");

    if idav <> fail then
        Print("  IdSmallGroup is available for this size. \n \n");
    else
        Print("  IdSmallGroup is not available for this size. \n \n");
    fi;
end );

#############################################################################
##
#F SMALL_GROUPS_INFORMATION[ 1 ]( size, smav, num )
##
SMALL_GROUPS_INFORMATION[ 1 ] := function( size, smav, num )
    Print("\n");
    Print("  The groups whose order factorises in at most 3 primes \n");
    Print("  have been classified by O. Hoelder. This classification is \n");
    Print("  used in the SmallGroups library. \n");
end;

SMALL_GROUPS_INFORMATION[ 2 ] := SMALL_GROUPS_INFORMATION[ 1 ];
SMALL_GROUPS_INFORMATION[ 3 ] := SMALL_GROUPS_INFORMATION[ 1 ];
SMALL_GROUPS_INFORMATION[ 4 ] := SMALL_GROUPS_INFORMATION[ 1 ];
SMALL_GROUPS_INFORMATION[ 5 ] := SMALL_GROUPS_INFORMATION[ 1 ];
SMALL_GROUPS_INFORMATION[ 6 ] := SMALL_GROUPS_INFORMATION[ 1 ];
SMALL_GROUPS_INFORMATION[ 7 ] := SMALL_GROUPS_INFORMATION[ 1 ];

#############################################################################
##
#F SMALL_GROUPS_INFORMATION[ 8 .. 10 ]( size, smav, num )
##
SMALL_GROUPS_INFORMATION[ 8 ] := function( size, smav, num )
    local ffid, prop, i, l;

    ffid := IdGroup( OneSmallGroup( size, FrattinifactorSize, size ) );
    prop := PROPERTIES_SMALL_GROUPS[ size ].frattFacs;

    if not IsPrimePowerInt( size ) then
        Print("  They are sorted by their Frattini factors. \n");
        i := 1;
        if ffid[ 2 ] > 1 then
            repeat
                if prop.pos[ i ][ 1 ] = -prop.pos[ i ][ 2 ] then
                    Print( "     ", prop.pos[ i ][ 1 ],
                       " has Frattini factor ", prop.frattFacs[ i ], ".\n"  );
                else
                    Print( "     ", prop.pos[ i ][ 1 ], " - ",
                       -prop.pos[ i ][ 2 ], " have Frattini factor ",
                       prop.frattFacs[ i ], ".\n"  );
                fi;
                i := i + 1;
            until prop.frattFacs[ i ] = ffid;
        fi;
        Print("     ", ffid[2], " - ", num,
              " have trivial Frattini subgroup.\n");
    else
        Print("  They are sorted by their ranks. \n");
        Print("     ", 1, " is cyclic. \n");
        i := 2;
        repeat
            l := Length( Factors( prop.frattFacs[ i ][1] ) );
            if prop.pos[ i ][ 1 ] = -prop.pos[ i ][ 2 ] then
                Print( "     ", prop.pos[ i ][ 1 ], " has rank ", l, ".\n"  );
            else
                Print( "     ", prop.pos[ i ][ 1 ], " - ",
                       -prop.pos[ i ][ 2 ], " have rank ", l, ".\n"  );
            fi;
            i := i + 1;
        until prop.frattFacs[ i ] = ffid;
        Print("     ", ffid[2], " is elementary abelian. \n");
    fi;

    Print( "\n  For the selection functions the values of the ",
           "following attributes \n  are precomputed and stored:\n ");
    if IsPrimePowerInt( size ) then
        Print( "    IsAbelian, PClassPGroup, RankPGroup,",
               " FrattinifactorSize and \n     FrattinifactorId. \n");
    else
        Print( "    IsAbelian, IsNilpotentGroup,",
               " IsSupersolvableGroup, IsSolvableGroup, \n     LGLength,",
               " FrattinifactorSize and FrattinifactorId. \n");
    fi;
end;
SMALL_GROUPS_INFORMATION[ 9 ] := SMALL_GROUPS_INFORMATION[ 8 ];
SMALL_GROUPS_INFORMATION[ 10 ] := SMALL_GROUPS_INFORMATION[ 8 ];

#############################################################################
##
#F SMALL_GROUPS_INFORMATION[ 11, 17 ]( size, smav, num )
##
SMALL_GROUPS_INFORMATION[ 11 ] := function( size, smav, num )
    local i, q;

    q := 2;
    if IsBound( smav.q ) then q := smav.q; fi;

    Print("  They are sorted by normal Sylow subgroups. \n");
    Print( "     1 - ", smav.pos[ 2 ], " are the nilpotent groups.\n" );
    for i in [ 2 .. Length( smav.types ) ] do
        Print( "     ", smav.pos[i] + 1, " - ", smav.pos[i+1] );
        if smav.types[ i ] = "p-autos" then
            Print( " have a normal Sylow ", q,"-subgroup. \n");
        elif smav.types[ i ] = "none-p-nil" then
            Print( " have no normal Sylow subgroup. \n");
        elif IsInt( smav.types[ i ] ) then
            Print( " have a normal Sylow ", smav.p, "-subgroup \n");
            Print( "                     with centralizer of index ");
            Print( q^smav.types[i],".\n");
        fi;
    od;
end;
SMALL_GROUPS_INFORMATION[ 17 ] := SMALL_GROUPS_INFORMATION[ 11 ];

#############################################################################
##
#F SMALL_GROUPS_INFORMATION[ 12 ]( size, smav, num )
##
SMALL_GROUPS_INFORMATION[ 12 ] := function( size, smav, num )

    if size = 1152 then
        Print("  They are sorted using Sylow subgroups. \n");
        Print("     1 - 2328 are nilpotent with Sylow 3-subgroup c9.\n" );
        Print("     2329 - 4656 are nilpotent with Sylow 3-subgroup 3^2.\n");
        Print("     4657 - 153312 are non-nilpotent with normal ");
        Print("Sylow 3-subgroup.\n");
        Print("     153313 - 157877 have no normal Sylow 3-subgroup.\n");
        return;
    fi;

    Print("  They are sorted using Hall subgroups. \n");
    Print( "     1 - 2328 are the nilpotent groups.\n" );
    Print( "     2329 - 236344 have a normal Hall (3,5)-subgroup.\n");
    Print( "     236345 - 240416 are solvable without normal Hall",
           " (3,5)-subgroup.\n");
    Print( "     240417 - 241004 are not solvable.\n" );
end;

#############################################################################
##
#F SMALL_GROUPS_INFORMATION[ 14 ]( size, smav, num )
##
SMALL_GROUPS_INFORMATION[ 14 ] := function( size, smav, num )

    Print( "     1 - 10494213 are the nilpotent groups.\n" );
    Print( "     10494214 - 408526597 have a normal Sylow 3-subgroup.\n" );
    Print( "     408526598 - 408544625 have a normal Sylow 2-subgroup.\n" );
    Print( "     408544626 - 408641062 have no normal Sylow subgroup.\n" );
end;

#############################################################################
##
#F SMALL_GROUPS_INFORMATION[ 18 ]( size, smav, num )
##
SMALL_GROUPS_INFORMATION[ 18 ] := function( size, smav, num )

    Print( "     1 is cyclic. \n");
    Print( "     2 - 10 have rank 2 and p-class 3.\n" );
    Print( "     11 - 386 have rank 2 and p-class 4.\n" );
    Print( "     387 - 1698 have rank 2 and p-class 5.\n" );
    Print( "     1699 - 2008 have rank 2 and p-class 6.\n" );
    Print( "     2009 - 2039 have rank 2 and p-class 7.\n" );
    Print( "     2040 - 2044 have rank 2 and p-class 8.\n" );
    Print( "     2045 has rank 3 and p-class 2.\n" );
    Print( "     2046 - 29398 have rank 3 and p-class 3.\n" );
    Print( "     29399 - 30617 have rank 3 and p-class 4.\n" );
    Print( "     30618 - 31239 have rank 3 and p-class 3.\n" );
    Print( "     31240 - 56685 have rank 3 and p-class 4.\n" );
    Print( "     56686 - 60615 have rank 3 and p-class 5.\n" );
    Print( "     60616 - 60894 have rank 3 and p-class 6.\n" );
    Print( "     60895 - 60903 have rank 3 and p-class 7.\n" );
    Print( "     60904 - 67612 have rank 4 and ", "p-class 2.\n" );
    Print( "     67613 - 387088 have rank 4 and ", "p-class 3.\n" );
    Print( "     387089 - 419734 have rank 4 and ", "p-class 4.\n" );
    Print( "     419735 - 420500 have rank 4 and ", "p-class 5.\n" );
    Print( "     420501 - 420514 have rank 4 and ", "p-class 6.\n" );
    Print( "     420515 - 6249623 have rank 5 and ", "p-class 2.\n" );
    Print( "     6249624 - 7529606 have rank 5 and ", "p-class 3.\n" );
    Print( "     7529607 - 7532374 have rank 5 and ", "p-class 4.\n" );
    Print( "     7532375 - 7532392 have rank 5 and ", "p-class 5.\n" );
    Print( "     7532393 - 10481221 have rank 6 and ", "p-class 2.\n" );
    Print( "     10481222 - 10493038 have rank 6 and ", "p-class 3.\n" );
    Print( "     10493039 - 10493061 have rank 6 and ", "p-class 4.\n" );
    Print( "     10493062 - 10494173 have rank 7 ", "and p-class 2.\n" );
    Print( "     10494174 - 10494200 have rank 7 ", "and p-class 3.\n" );
    Print( "     10494201 - 10494212 have rank 8 ", "and p-class 2.\n" );
    Print( "     10494213 is elementary abelian.\n");
end;

#############################################################################
##
#F SMALL_GROUPS_INFORMATION[ 19 ]( size, smav, num )
##
SMALL_GROUPS_INFORMATION[ 19 ] := function( size, smav, num )

  Print("  They are sorted by their ranks. \n");
  Print( "     1 is cyclic. \n");
  Print( "     2 - 10 have rank 2. \n");
  Print( "     11 - 14 have rank 3. \n");
  Print( "     15 is elementary abelian. \n");
end;

#############################################################################
##
#F SMALL_GROUPS_INFORMATION[ 20 ]( size, smav, num )
##
SMALL_GROUPS_INFORMATION[ 20 ] := function( size, smav, num )
    local p, a, b, c;

    p := Factors(size)[1];
    a:=27 + p   + 2*GcdInt(p-1,3) + GcdInt(p-1,4);
    b:=54 + 2*p + 2*GcdInt(p-1,3) + GcdInt(p-1,4);
    c:=60 + 2*p + 2*GcdInt(p-1,3) + GcdInt(p-1,4);

    Print( "  They are sorted by their ranks.\n" );
    Print( "     1 is cyclic.\n");
    Print( "     2 - ",a," have rank 2. \n");
    Print( "     ",a+1," - ",b," have rank 3. \n");
    Print( "     ",b+1," - ",c," have rank 4. \n");
    Print( "     ",c+1," is elementary abelian. \n");
end;

#############################################################################
##
#F SMALL_GROUPS_INFORMATION[ 21 ]( size, smav, num )
##
SMALL_GROUPS_INFORMATION[ 21 ] := function( size, smav, num )

   Print( " \n");
   Print( "      Easterfield (1940) constructed a list of the groups of\n");
   Print( "      order p^6 for p >= 5.\n \n");

   Print( "      The database of parametrised presentations for the groups \n");
   Print( "      with order p^6 for p >= 5 is based on the Easterfield \n");
   Print( "      list, corrected by Newman, O'Brien and Vaughan-Lee (2004).\n");
   Print( "      It differs only in the addition of groups in isoclinism \n");
   Print( "      family $\\Phi_{13}$, in using the James (1980) presentations \n");
   Print( "      for the groups in $\\Phi_{19}$, and a small number of \n");
   Print( "      typographical amendments. The linear ordering employed is \n");
   Print( "      very close to that of Easterfield. \n \n");

   Print( "      Each group with order $p^6$ is described by a power- \n");
   Print( "      commutator presentation on 6 generators and 21 relations:\n");
   Print( "      15 are commutator relations and 6 are power relations. \n");
   Print( "      Each presentation has the prime $p$ as a parameter. \n");
   Print( "      The database contains about 500 parametrised \n");
   Print( "      presentations, most of these have $p$ as the only \n");
   Print( "      parameter. \n");

end;

#############################################################################
##
#F SMALL_GROUPS_INFORMATION[ 24 ]( size, smav, num )
##
SMALL_GROUPS_INFORMATION[ 24 ] := function( size, smav, num )
    local i, set, c;

    Print( "\n" );
    Print( "  The groups of squarefree order have a cylic socle and a " );
    Print( "cylic socle factor.\n" );

    Print( "\n" );
    i := 0;
    for set in smav.sets do
        c := Product( smav.primes{ set.kp } );
        if c = 1 then
            Print( "    1 is abelian\n" );
        elif set.number = 1 then
            Print( "    ", i + 1, " has socle C_" );
            Print( size / c, " and factor C_", c, "\n" );
        else
            Print( "    ", i + 1, " - ", i + set.number, " have socle C_" );
            Print( size / c, " and factor C_", c, "\n" );
        fi;
        i := i + set.number;
    od;
end;

#############################################################################
##
#F SMALL_GROUPS_INFORMATION[ 25 ]( size, smav, num )
##
SMALL_GROUPS_INFORMATION[ 25 ] := function( size, smav, num )
    local i, set, c;
    Print( "\n" );
    Print( "  The groups of cubefree order are either solvable or a direct ",
           "product of \n  the form PSL( 2, p ) x solvable group. ",
           "The cubefree solvable groups are \n  determined by their Frattini",
           " factor.\n\n" );

    i := 0;
    for set in smav.sets do
      if set.psl_p = 1 then
        if set.size_phi = 1 then
          if set.number = 1 then
            Print( "    ", i + 1, " is solvable and Frattini free\n" );
          else
            Print( "    ", i + 1, " - ", i + set.number, " are solvable ",
                   "and Frattini free\n" );
          fi;
        else
          if set.number = 1 then
            Print( "    ", i + 1, " is solvable with Frattini factor of ",
                   "size ", set.size_ff, "\n" );
          else
            Print( "    ", i + 1, " - ", i + set.number, " are solvable ",
                   "with Frattini factor of size ", set.size_ff, "\n" );
          fi;
        fi;
      elif
        set.size_ff = 1 then
          Print( "    ", i + 1, " is PSL( 2, ", set.psl_p, " )\n" );
      else
        if set.size_phi = 1 then
          if set.number = 1 then
            Print( "    ", i + 1, " is PSL( 2, ", set.psl_p, " ) x F, F ",
                   "solvable and Frattini free of order ", set.size_ff, "\n");
          else
            Print( "    ", i + 1, " - ", i + set.number, " are PSL( 2, ",
                   set.psl_p, " ) x F_i, F_i solvable ",
                   "Frattini free of order ", set.size_ff, "\n" );
          fi;
        else
          if set.number = 1 then
            Print( "    ", i + 1, " is PSL( 2, ", set.psl_p, " ) x G, G ",
                   "solvable of order ", set.size_ff * set.size_phi,
                   " with a Frattini factor\n      of order ", set.size_ff,
                   "\n");
          else
            Print( "    ", i + 1, " - ", i + set.number, " are PSL( 2, ",
                   set.psl_p, " ) x G_i, G_i ", "solvable of order ",
                   set.size_ff * set.size_phi, " with a",
                   "\n      Frattini factor of order ", set.size_ff, "\n");
          fi;
        fi;
      fi;
      i := i + set.number;
    od;
end;

#############################################################################
##
#F SMALL_GROUPS_INFORMATION[ 26 ]( size, smav, num )
##
SMALL_GROUPS_INFORMATION[ 26 ] := function( size, smav, num )

   Print( " \n");
   Print( "      E.A. O'Brien and M.R. Vaughan-Lee determined presentations\n");
   Print( "      of the groups with order p^7. A preprint of their paper is\n");
   Print( "      available at\n" );
   Print( "      http://www.math.auckland.ac.nz/%7Eobrien/research/p7/paper-p7.pdf\n\n" );
   Print( "      For p in { 3, 5, 7, 11 } explicit lists of groups of order\n");
   Print( "      p^7 have been produced and stored into the database.\n\n");
   Print( "      Giving the power commutator presentations of any of these\n");
   Print( "      groups using a standard notation they might be reduced to 35\n");
   Print( "      elements of the group or a 245 p-digit number.\n\n");
   Print( "      Only 56 of these digits may be unlike 0 for any group and\n");
   Print( "      even these 56 digits are mostly like 0. Further on these\n");
   Print( "      digits are often quite likely for sequences of subsequent\n");
   Print( "      groups. Thus storage of groups was done by finding a so\n");
   Print( "      called head group and a so called tail. Along the tail\n");
   Print( "      only the different digits compared to the head are relevant.\n");
   Print( "      Even the tails occur more or less often and this is used\n");
   Print( "      to improve storage too. Since p^7 is too big the data is\n");
   Print( "      stored into some remaing holes of SMALL_GROUP_LIB at\n");
   Print( "      Primes[ p + 10 ].\n");

end;

max_stars_repo_name: hulpke/smallgrp

max_issues_repo_name: hulpke/smallgrp

max_forks_repo_name: gap-packages/smallgrp


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

Class(CellCodegen, VectorCodegen, rec(
    Formula := meth(self, o, y, x, opts)
        local icode, datas, prog, params, sub, initsub, io;
        if IsBound(opts.XType) and not IsList(x) then
            x.t := TPtr(opts.XType);
            if IsBound(opts.useRestrict) and opts.useRestrict then
                x.t := x.t.restrict();
            fi;
        fi;
        if IsBound(opts.YType) and not IsList(y) then
            y.t := TPtr(opts.YType);
            if IsBound(opts.useRestrict) and opts.useRestrict then
                y.t := y.t.restrict();
            fi;
        fi;

        o := o.child(1);
        params := Set(Collect(o, param));

        datas := Collect(o, FDataOfs);
        o := BlockSums(opts.globalUnrolling, o);
        icode := ESReduce(self(o, y, x, opts),opts);
        icode := RemoveAssignAcc(icode);
        #Error("BP");
        icode := BlockUnroll(icode, opts);
        # icode := PowerOpt(icode);
        icode := DeclareHidden(icode);
        if IsBound(opts.isFixedPoint) and opts.isFixedPoint then
            icode := FixedPointCode(icode, opts.bits, opts.fracbits);
        fi;

        # Insert sw pipelining here
        icode := MarkPreds(icode);
        icode := MarkDefUse(icode);
        #SubstTopDown(icode, loop, e->MarkSWPLoops(e));
        #SubstTopDown(icode, loop_sw, e->SoftwarePipeline(e));

        io := When(x=y, [x], [y, x]);
        sub := Cond(IsBound(opts.subName), opts.subName, "transform");
        initsub := Cond(IsBound(opts.subName), Concat("init_", opts.subName), "init");
        icode := func(TVoid, sub, Concatenation(io, params), icode);

        if IsBound(opts.generateInitFunc) and opts.generateInitFunc then
            prog := program(
                decl(List(datas, x->x.var),
                    chain(
                        func(TVoid, initsub, params, chain(List(datas, x -> SReduce(x.var.init, opts)))),
                        icode
                    )));
        else
            prog := program( func(TVoid, initsub, params, chain()), icode);
        fi;
        prog.dimensions := o.dims();
        return prog;
    end,

));

ext: gi


text:
%Terminals
    DollarSign ::= '$'
    Percent ::= '%'
    _
    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
%End

%Headers
    /.
        //
        // Each upper case letter is mapped into its corresponding
        // lower case counterpart. For example, if an 'A' appears
        // in the input, it is mapped into $sym_type::$prefix$a$suffix$ just
        // like 'a'.
        //
         inline static int tokenKind[128] = {};

        static bool static_init()
        {
            tokenKind['$'] = $sym_type::$prefix$DollarSign$suffix$;
            tokenKind['%'] = $sym_type::$prefix$Percent$suffix$;
            tokenKind['_'] = $sym_type::$prefix$_$suffix$;

            tokenKind['a'] = $sym_type::$prefix$a$suffix$;
            tokenKind['b'] = $sym_type::$prefix$b$suffix$;
            tokenKind['c'] = $sym_type::$prefix$c$suffix$;
            tokenKind['d'] = $sym_type::$prefix$d$suffix$;
            tokenKind['e'] = $sym_type::$prefix$e$suffix$;
            tokenKind['f'] = $sym_type::$prefix$f$suffix$;
            tokenKind['g'] = $sym_type::$prefix$g$suffix$;
            tokenKind['h'] = $sym_type::$prefix$h$suffix$;
            tokenKind['i'] = $sym_type::$prefix$i$suffix$;
            tokenKind['j'] = $sym_type::$prefix$j$suffix$;
            tokenKind['k'] = $sym_type::$prefix$k$suffix$;
            tokenKind['l'] = $sym_type::$prefix$l$suffix$;
            tokenKind['m'] = $sym_type::$prefix$m$suffix$;
            tokenKind['n'] = $sym_type::$prefix$n$suffix$;
            tokenKind['o'] = $sym_type::$prefix$o$suffix$;
            tokenKind['p'] = $sym_type::$prefix$p$suffix$;
            tokenKind['q'] = $sym_type::$prefix$q$suffix$;
            tokenKind['r'] = $sym_type::$prefix$r$suffix$;
            tokenKind['s'] = $sym_type::$prefix$s$suffix$;
            tokenKind['t'] = $sym_type::$prefix$t$suffix$;
            tokenKind['u'] = $sym_type::$prefix$u$suffix$;
            tokenKind['v'] = $sym_type::$prefix$v$suffix$;
            tokenKind['w'] = $sym_type::$prefix$w$suffix$;
            tokenKind['x'] = $sym_type::$prefix$x$suffix$;
            tokenKind['y'] = $sym_type::$prefix$y$suffix$;
            tokenKind['z'] = $sym_type::$prefix$z$suffix$;

            tokenKind['A'] = $sym_type::$prefix$a$suffix$;
            tokenKind['B'] = $sym_type::$prefix$b$suffix$;
            tokenKind['C'] = $sym_type::$prefix$c$suffix$;
            tokenKind['D'] = $sym_type::$prefix$d$suffix$;
            tokenKind['E'] = $sym_type::$prefix$e$suffix$;
            tokenKind['F'] = $sym_type::$prefix$f$suffix$;
            tokenKind['G'] = $sym_type::$prefix$g$suffix$;
            tokenKind['H'] = $sym_type::$prefix$h$suffix$;
            tokenKind['I'] = $sym_type::$prefix$i$suffix$;
            tokenKind['J'] = $sym_type::$prefix$j$suffix$;
            tokenKind['K'] = $sym_type::$prefix$k$suffix$;
            tokenKind['L'] = $sym_type::$prefix$l$suffix$;
            tokenKind['M'] = $sym_type::$prefix$m$suffix$;
            tokenKind['N'] = $sym_type::$prefix$n$suffix$;
            tokenKind['O'] = $sym_type::$prefix$o$suffix$;
            tokenKind['P'] = $sym_type::$prefix$p$suffix$;
            tokenKind['Q'] = $sym_type::$prefix$q$suffix$;
            tokenKind['R'] = $sym_type::$prefix$r$suffix$;
            tokenKind['S'] = $sym_type::$prefix$s$suffix$;
            tokenKind['T'] = $sym_type::$prefix$t$suffix$;
            tokenKind['U'] = $sym_type::$prefix$u$suffix$;
            tokenKind['V'] = $sym_type::$prefix$v$suffix$;
            tokenKind['W'] = $sym_type::$prefix$w$suffix$;
            tokenKind['X'] = $sym_type::$prefix$x$suffix$;
            tokenKind['Y'] = $sym_type::$prefix$y$suffix$;
            tokenKind['Z'] = $sym_type::$prefix$z$suffix$;
            return true;
        };
        inline static bool ddddd = static_init();
         int getKind(wchar_t c)
        {
            return (c < 128 ? tokenKind[c] : 0);
        }
    ./
%End

max_stars_repo_name: kuafuwang/LPG-VScode

max_issues_repo_name: kuafuwang/LPG-VScode

max_forks_repo_name: kuafuwang/LPG-VScode


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

Class(ScriptGenSSE, ScriptGenBase, rec(

	_arch := () -> "SSE",

	_init := meth(self)
		self._setTransform(SGKEY_FFT);
		self._setType(SGKEY_SPCX);
		self._setSize(16);
		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
			szs := Filtered([1..512], i->ForAll(Factors(i), j->j<=16) and IsInt(i/4));
			Append(szs, Filtered([513..1024], i->ForAll(Factors(i), j->j<=16) and IsInt(i/16)));
			Append(szs, List([11..16], i->2^i));
		elif xform in [SGKEY_FFT_2D, SGKEY_IFFT_2D] then
			szs := List(Filtered([1..360], i->ForAll(Factors(i), j->j<=16) and IsInt(i/4)), n -> [n, n]);
		elif xform = SGKEY_WHT then
			if type = SGKEY_DPCX then
				szs := List([2..10], i->2^i);
			else
				szs := List([4..10], 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;

		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    := "SSE_4x32f";
			use_cx := ForAny(szs, s -> not IsInt(s / 16));
		else
			isa    := "SSE_2x64f";
			use_cx := false;
		fi;

		if xform in [SGKEY_FFT, SGKEY_IFFT] then
			optrec := rec (
				tsplRader:=false,
				tsplBluestein:=false,
				tsplPFA:=false,
				oddSizes:=false,
				interleavedComplex := true,
			);
			if use_cx then
				if (type = SGKEY_SPCX) then
					optrec.cplxVect := true;
					optrec.realVect := false;
					optrec.PRDFT := false;
					optrec.URDFT := 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 (
				svct := true,
				stdTTensor := true,
				tsplPFA := false
			);
			if (type = SGKEY_SPCX) then
				optrec.oddSizes := false;
			else
				optrec.verify := true;
				optrec.oddSizes := true;
			fi;
			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,

)); # Class ScriptGenSSE

SetScriptGenConstructor(ScriptGenSSE);


text:
MersenneSmallFactor := function(n)
    local k, m, d;
    if IsPrime(n) then
        d := 2*n;
        m := 1;
        for k in [1 .. 1000000] do
            m := m + d;
            if PowerModInt(2, n, m) = 1 then
                return m;
            fi;
        od;
    fi;
    return fail;
end;

# If n is not prime, fail immediately
MersenneSmallFactor(15);
# fail

MersenneSmallFactor(929);
# 13007

MersenneSmallFactor(1009);
# 3454817

# We stop at k = 1000000 in 2*k*n + 1, so it may fail if 2^n - 1 has only larger factors
MersenneSmallFactor(101);
# fail

FactorsInt(2^101-1);
# [ 7432339208719, 341117531003194129 ]

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

# This file implements a register allocator for straight line code
# Based on Belady's farthest-first eviction strategy.  For normal
# three-operand code without memory ops, this should yield the optimal
# (w.r.t. number of spill-related load-stores) code.
#
# For two-operand code without memory ops, this can still be optimal,
# but I'm not sure.
#
# For two-operand code *with* memory ops (=x86) it is not clear at
# all, whether this simple allocator is best. Main problem is that we
# have a choice of spilling versus using memop in an instruction.
#
# The output of this pass should be fed into C-compiler, and spill
# space here is represented as a regular array. The hope is that
# C-compiler won't touch the spill array (S), and therefore will not
# run out of registers.
#
# Potential problems: on 2-operand architectures (x86), the compiler
# can internally operate on 3-operand code, and then convert to
# 2-operand at some point. This may render this allocator useless,
# because 3->2 conversion is likely to destroy lot's of information.
# NOTE: Investigate this further.
#

# we are assuming SSA code
MarkLiveness := function(code)
    local live, add, remove, v, c, already_dead, n;
    MarkDefUse(code);

    # Starting from the statement, invariant : variable is dead afters its last use
    already_dead := Set([]);
    for c in Reversed(code.cmds) do
        c.dead := Difference(PredCmd(c), already_dead);
	UniteSet(already_dead, c.dead);
    od;

    live := Set([]);
    n := 1;
    for c in code.cmds do
        Constraint(IsAssign(c));
	if SuccCmd(c) <> [] then AddSet(live, c.loc); fi;
	SubtractSet(live, c.dead);
	c.live := Copy(live);
	c.n := n; # position in a chain
	n := n+1;
    od;
    return code;
end;

# RequiredVars(<chain>) - minimum number of required variables
# MarkLiveness should be run beforehand
#
RequiredVars := (code) -> Maximum0(List(code.cmds, c->Length(c.live)));

# ReuseVars(<chain>) - reusing variables in the chain
# MarkLiveness should be run beforehand,
# output is no longer in SSA form
#

ReuseVars := function(code)
    local free, current, map, reduced, new, out, v, c, f, alive;
    reduced := [];
    free    := [];
    current := [];

    map     := tab();

    for c in code.cmds do

        alive := Set(c.live :: Filtered(c.dead, e -> IsVar(e) and not IsArrayT(e.t) and not IsPtrT(e.t)));
        new   := Difference(alive, current);
        out   := Difference(current, alive);

        for v in out do
            if IsBound(map.(v.id)) then
                free := [map.(v.id)] :: free;
                Unbind(map.(v.id));
            else
                free := free :: [v];
            fi;
        od;
        if Length(new)=1 then
            f := PositionProperty(free, e -> e.t=new[1].t);
            if f<>false then
                map.(new[1].id) := free[f];
                free := ListWithout(free, f);
            fi;
        fi;

        Add(reduced, SubstParamsCustom(c, map, [var]));
        current := alive;
    od;

    return chain(reduced);
end;

# RequiredRegs(<chain>) - minimum number of required registers
# MarkLiveness should be run beforehand
#
RequiredRegs := (code, t) -> Maximum(List(code.cmds, c->Length(
	    Filtered(c.live, x->IsBound(x.t) and x.t=t))));

# spill using Belady's MIN heuristic - spill the value used farthest away
_spill := function(n, regmap, regvals, S, free_slots, newcmds)
    local dists, maxdist, spill, evict, r;

    dists := List(regvals, v -> let(
	    dd := Filtered(List(SuccLoc(v), succ -> succ.def.n), succn -> succn >= n),
            # if successor list (beyond current n) is empty, the variable is dead and should
	    # be evicted, assign it a really large dist value 2^20
	    When(dd=[], 2^20, Minimum(dd))));

    maxdist := Maximum(dists);
    evict := regvals[Position(dists, maxdist)];
    r := regmap.(evict.id); # freed up register

    if not IsBound(evict.value) then # do not save constants into spill slot
	Constraint(Length(free_slots) > 0);
	spill := nth(S, free_slots[1]);
	RemoveSet(free_slots, free_slots[1]);
	Add(newcmds, assign(spill, r));
	regmap.(evict.id) := spill;
    else
	Unbind(regmap.(evict.id));
    fi;

    RemoveSet(regvals, evict);

    return r;
end;

_getreg := function(n, newvar, free, regmap, regvals, S, free_slots, newcmds)
    local r;
    if Length(free) > 0 then
        # we have free registers, grab <r> from free list
	r := free[1];
	RemoveSet(free, r);
    else
 	# no free registers spill something to free up <r>
	r := _spill(n, regmap, regvals, S, free_slots, newcmds);
    fi;
    AddSet(regvals, newvar);
    regmap.(newvar.id) := r;
    return r;
end;

_reload := function(n, vars, free, regmap, regvals, S, free_slots, newcmds)
    local v, r, reload_from;
    for v in vars do
        reload_from := When(IsBound(regmap.(v.id)), regmap.(v.id), v);
        r := _getreg(n, v, free, regmap, regvals, S, free_slots, newcmds);
	Add(newcmds, assign(r, reload_from));
        # if <reload_from> is a spill location (ie not a variable),
        # then return the  spill slot to the free list
	if not IsVar(reload_from) then
	AddSet(free_slots, reload_from.idx); fi;
    od;
end;

# MarkLiveness should be run beforehand
#
RegAlloc := function(code, numregs, t)
    local regs, stack, regmap, regvals, c, cc, free, dead, r, v,
         dists, maxdist, evict, newcmds,
	 S, free_slots, reqregs, spill, deadregs,
	 reloads;

    reqregs := RequiredRegs(code, t);
    S := var("S", TArray(t, When(reqregs-numregs+1 > 0, reqregs-numregs+21, 0)));
    # free spill slots
    free_slots := Set([0..S.t.size-1]);

    regmap := rec();
    free := Set(List([0..numregs-1], x -> var(Concat("r", StringInt(x)), t)));
    for r in free do r.isReg := true; od;

    regvals := Set([]); # set of all variables that reside in registers (=not spilled)

    newcmds := [];
    for c in code.cmds do
        Constraint(IsAssign(c));
	# plug in allocated registers instead of SSA vars
	cc := Copy(c);

	reloads := Union(
	    # do not reload constants (in register for which regmap.(..) == var, and external)
	    Filtered(cc.op_inout(), x->IsVar(x) and (not IsBound(regmap.(x.id)) or # external
		                                     not IsVar(regmap.(x.id))) and not (TPtr = ObjId(x.t) or TInt = x.t)), # or spilled

	    # do not reload constants or other inputs used only once
	    # ideally we should use Belady heuristic again
	    Filtered(cc.op_in(), x->IsVar(x) and (not IsBound(regmap.(x.id)) or
		                                  not IsVar(regmap.(x.id))
	                                     and (not Length(SuccCmd(x.def))=1)) and not (TPtr = ObjId(x.t) or  TInt = x.t)));

	_reload(cc.n, reloads, free, regmap, regvals, S, free_slots, newcmds);
	if (cc.loc in cc.op_inout()) or (cc.loc in cc.op_in()) then
	    cc := SubstVars(cc, regmap);
	else
	    cc.exp := SubstVars(cc.exp, regmap);
	fi;
 	# allocate a register if c.loc only is a variable of type t
	if IsVar(cc.loc) and not IsBound(cc.loc.isReg) and cc.loc.t = t then
	    deadregs := Filtered(cc.dead, v -> IsVar(v) and v in regvals);
	    if deadregs <> [] then
		# grab reg from the dead list (this automatically implements coalescing)
		v := deadregs[1];

		RemoveSet(regvals, v);
		AddSet(regvals, cc.loc);
		r := regmap.(v.id);
		Unbind(regmap.(v.id));
		regmap.(cc.loc.id) := r;
		cc.loc := r;
	    else
		r := _getreg(cc.n, cc.loc, free, regmap, regvals, S, free_slots, newcmds);
		cc.loc := r;
	    fi;
	fi;
	# do not output dummy moves (ie assign(r1, r1)).
	if not (ObjId(cc)=assign and cc.loc = cc.exp) then
	    Add(newcmds, cc);
	fi;

	# put registers associated with dead variables on the free list
	for v in cc.dead do
	    if IsVar(v) and (v in regvals) then
		AddSet(free, regmap.(v.id));
		RemoveSet(regvals, v);
		Unbind(regmap.(v.id));
	    fi;
	od;

    od;
    return chain(newcmds);
end;


text:
# Creating an array of 0
a := NullMat(2, 2);
# [ [ 0, 0 ], [ 0, 0 ] ]

# Some assignments
a[1][1] := 4;
a[1][2] := 5;
a[2][1] := 3;
a[2][2] := 4;

a
# [ [ 4, 5 ], [ 3, 4 ] ]

Determinant(a);
# 1


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

Class(SKLR_Vx1i, SIMD_ISA, rec(

    includes     := () -> ["<stdlib.h>"] :: _MM_MALLOC(),
    active       := true,
    ctype        := "char",
    instr        := [],
    bits         := 1,
    isFloat      := false,
    isFixedPoint := false,
    splopts      := rec(),
    alignment    := 8,

    autolib := rec(
        includes := () -> ["<include/sp_bits.h>"],
        includesTimer := () -> [],
    ),

    _op_load_u  := abstract(),
    _op_bcast   := abstract(),
    _op_store_u := abstract(),

    dupload  := (self, y, x) >> Checked(ObjId(x)=nth, # NOTE: is there a better way?
	let(base := x.loc,
	    ofs  := x.idx,
	    v    := self.v,
	    xvec := Cond(IsUnalignedPtrT(base.t), self._op_load_u(base, idiv(ofs,v)*v, v),
		         vtref(self.t, base, idiv(ofs, v))),
	    assign(y, self._op_bcast(xvec, imod(ofs, v))))),

    svload := [ [ ], # load using subvecs of len 1
                [ ], # load using subvecs of len 2
                [ ], # load using subvecs of len 4
    ],

    svstore := [ [ ], # store using subvecs of len 1
                 [ ], # store using subvecs of len 2
                 [ ], # store using subvecs of len 4
    ],

    # keep the n lower scalars and zero the other ones
    mask_l := (self, c, n) >> Cond( n = self.v, c,
        bin_and(c, self.val(Replicate(n, 1) :: Replicate(self.v - n, 0)))),
    mask_h := (self, c, n) >> Cond( n = self.v, c,
        bin_and(c, self.val(Replicate(n, 0) :: Replicate(self.v - n, 1)))),

    loadCont := (self, n, y, yofs, x, xofs, xofs_align, opts) >> let(
	a  := _unwrap(xofs_align),
	nn := _unwrap(n),
	yy := vtref(self.t, y, yofs),
	m  := x -> self.mask_l(x, nn),
	Cond(a = 0 and not IsUnalignedPtrT(x.t),
	         assign(yy, m(vtref(self.t, x, xofs/self.v))),

	     # known alignment, sv is small, so that we only need 1 aligned load + 1 shift + mask
	     ((IsInt(a) and (nn <= self.v - a)) or nn=1) and not IsUnalignedPtrT(x.t),
		 let(v1 := vtref(self.t, x, idiv(xofs, self.v)),
		     assign(yy, m(bin_shr(v1, a)))),

	     # known alignment, sv covers 2 vectors, 2 aligned loads + 2 shifts + mask
	     # NB: no masking is needed because shifts will do the job
	     IsInt(a) and not IsUnalignedPtrT(x.t),
		 let(v1 := vtref(self.t, x, (xofs - a)/self.v),
		     v2 := vtref(self.t, x, (xofs - a)/self.v + 1),
		     assign(yy, m(bin_or(bin_shr(v1, a), bin_shl(v2, self.v - a))))),
             # else, unknown alignment, use unaligned load
             assign(yy, self._op_load_u(x, xofs, nn))
        )
    ),

    storeCont := (self, n, y, yofs, yofs_align, x, xofs, opts) >> let(
	a  := _unwrap(yofs_align),
	nn := _unwrap(n),
	xx := vtref(self.t, x, xofs),
	yy := vtref(self.t, y, yofs/self.v),
	Cond(nn = self.v and a = 0 and not IsUnalignedPtrT(y.t),
	        assign(yy, xx),
	     #else
	        self._op_store_u(y, yofs, xx, nn))),

    rotate_left := (self, shift) >> ((y, x) -> assign(vtref(self.t, y, 0), rCyclicShift(vtref(self.t, x, 0), shift, self.v))),

    kswap := (self, y, x, k, mask) >> let( u := var.fresh_t("U", self.t),
                                         chain(assign( u, bin_and(bin_xor(x, bin_shr(x, 2^(k-1))), self.t.value(mask))),
                                               assign( y, bin_xor(bin_xor(x, u), bin_shl(u, 2^(k-1)))))),

    kexch := (self, y1, y2, x1, x2, mask) >> let( u := var.fresh_t("U", self.t),
                                         chain(assign(  u, bin_and(bin_xor(x1, x2), self.t.value(mask))),
                                               assign( y1, bin_xor(x1, u)),
                                               assign( y2, bin_xor(x2, u)))),
));

Class(SKLR_16x1i, SKLR_Vx1i, rec(
    # countrec below is invalid
    countrec := rec(
        ops := [
            [add, sub],
	    [mul],
            [sklr_bcast_16x1i], # shuffles
            [sklr_loadu_16x1i, sklr_storeu_16x1i],
            [deref],
            Value      # Value without [] is a keyword in countOps !!
        ],
        printstrings := ["[adds]", "[mults]", "[vperms]", "[svldst]", "[vldst]", "[vval]"],
        type := "TVect",
        arithcost := (self, opcount) >> opcount[1]+opcount[2]
    ),

    info     := "Scalar 16 x 1-bit",
    v        := 16,
    t        := BitVector(16),

    v_ones   := BitVector(16).one(),
    v_zeros  := BitVector(16).zero(),
    val      := bits -> BitVector(16).value(bits),

    _op_load_u  := (self, ptr, offs, elts)      >> sklr_loadu_16x1i(ptr, offs, elts),
    _op_bcast   := (self, loc, elt_num)         >> sklr_bcast_16x1i(loc, elt_num),
    _op_store_u := (self, ptr, offs, src, elts) >> sklr_storeu_16x1i(ptr, offs, src, elts),
));

Class(SKLR_32x1i, SKLR_Vx1i, rec(
    # countrec below is invalid
    countrec := rec(
        ops := [
            [add, sub],
	    [mul],
            [sklr_bcast_32x1i], # shuffles
            [sklr_loadu_32x1i, sklr_storeu_32x1i],
            [deref],
            Value      # Value without [] is a keyword in countOps !!
        ],
        printstrings := ["[adds]", "[mults]", "[vperms]", "[svldst]", "[vldst]", "[vval]"],
        type := "TVect",
        arithcost := (self, opcount) >> opcount[1]+opcount[2]
    ),

    info     := "Scalar 32 x 1-bit",
    v        := 32,
    t        := BitVector(32),

    v_ones   := BitVector(32).one(),
    v_zeros  := BitVector(32).zero(),
    val      := bits -> BitVector(32).value(bits),

    _op_load_u  := (self, ptr, offs, elts)      >> sklr_loadu_32x1i(ptr, offs, elts),
    _op_bcast   := (self, loc, elt_num)         >> sklr_bcast_32x1i(loc, elt_num),
    _op_store_u := (self, ptr, offs, src, elts) >> sklr_storeu_32x1i(ptr, offs, src, elts),

));

Class(SKLR_64x1i, SKLR_Vx1i, rec(
    # countrec below is invalid
    countrec := rec(
        ops := [
            [add, sub],
	    [mul],
            [sklr_bcast_64x1i], # shuffles
            [sklr_loadu_64x1i, sklr_storeu_64x1i],
            [deref],
            Value      # Value without [] is a keyword in countOps !!
        ],
        printstrings := ["[adds]", "[mults]", "[vperms]", "[svldst]", "[vldst]", "[vval]"],
        type := "TVect",
        arithcost := (self, opcount) >> opcount[1]+opcount[2]
    ),

    info     := "Scalar 64 x 1-bit",
    v        := 64,
    t        := BitVector(64),

    v_ones   := BitVector(64).one(),
    v_zeros  := BitVector(64).zero(),
    val      := bits -> BitVector(64).value(bits),

    _op_load_u  := (self, ptr, offs, elts)      >> sklr_loadu_64x1i(ptr, offs, elts),
    _op_bcast   := (self, loc, elt_num)         >> sklr_bcast_64x1i(loc, elt_num),
    _op_store_u := (self, ptr, offs, src, elts) >> sklr_storeu_64x1i(ptr, offs, src, elts),
));

RewriteRules(RulesStrengthReduce, rec(
    aligned_loadu_32x1 := Rule( @(1, sklr_loadu_32x1i, x -> IsInt(_unwrap(x.args[2] mod 32)) and not IsUnalignedPtrT(x.args[1])),
                             e -> let(
                                     offs := imod(e.args[2], 32),
                                     xx0  := vtref(e.t, e.args[1], idiv(e.args[2], 32)),
                                     xx1  := vtref(e.t, e.args[1], idiv(e.args[2], 32)+1),
                                     nn   := _unwrap(e.args[3]),
                                     bin_and( When(offs + nn <= 32,
                                             bin_shr(xx0, offs),
                                             bin_or(bin_shr(xx0, offs), bin_shl(xx1, 32 - offs))),
                                         e.t.value(Replicate(nn, 1) :: Replicate(32 - nn, 0)))
                                     )),
    aligned_loadu_64x1 := Rule( @(1, sklr_loadu_64x1i, x -> IsInt(_unwrap(x.args[2] mod 64)) and not IsUnalignedPtrT(x.args[1])),
                             e -> let(
                                     offs := imod(e.args[2], 64),
                                     xx0  := vtref(e.t, e.args[1], idiv(e.args[2], 64)),
                                     xx1  := vtref(e.t, e.args[1], idiv(e.args[2], 64)+1),
                                     nn   := _unwrap(e.args[3]),
                                     bin_and( When(offs + nn <= 64,
                                             bin_shr(xx0, offs),
                                             bin_or(bin_shr(xx0, offs), bin_shl(xx1, 64 - offs))),
                                         e.t.value(Replicate(nn, 1) :: Replicate(64 - nn, 0)))
                                     )),

));

Class(SKLR_32x1i_to_SSE_16x8i, ISA_Bridge, rec(
    isa_from    := SKLR_32x1i,
    isa_to      := SSE_16x8i(T_Int(8)),

    code        := (self, y, x, opts) >> let(
                    xx := (offs) -> vtref(self.isa_from.t, x, offs),
                    yy := (offs) -> vtref(self.isa_to.t,   y, offs),
                    a  := var.fresh_t("U", T_UInt(32)),
                    b0 := var.fresh_t("U", T_UInt(32)),
                    b1 := var.fresh_t("U", T_UInt(32)),
                    b2 := var.fresh_t("U", T_UInt(32)),
                    b3 := var.fresh_t("U", T_UInt(32)),
                    mask := T_UInt(32).value(1 + 256 + 65536 + 16777216),
                    decl([a,b0,b1,b2,b3], chain(
                        assign( a, tcast(a.t, xx(0)) ),
                        assign( b0, bin_and(            a, mask)),
                        assign( b1, bin_and(bin_shr(a, 1), mask)),
                        assign( b2, bin_and(bin_shr(a, 2), mask)),
                        assign( b3, bin_and(bin_shr(a, 3), mask)),
                        assign( yy(0), tcast(self.isa_to.t, vpack(b0, b1, b2, b3))),
                        assign( b0, bin_and(bin_shr(a, 4), mask)),
                        assign( b1, bin_and(bin_shr(a, 5), mask)),
                        assign( b2, bin_and(bin_shr(a, 6), mask)),
                        assign( b3, bin_and(bin_shr(a, 7), mask)),
                        assign( yy(1), tcast(self.isa_to.t, vpack(b0, b1, b2, b3)))
                    ))),

    toAMat := self >> L(self.isa_from.v, 8).toAMat(),
    toSpl  := self >> Cvt(self)*TL(self.isa_from.v, div(self.isa_from.v, 8)).withTags([AVecReg(self.isa_from)])
));

Class(SKLR_64x1i_to_SSE_16x8i, SKLR_32x1i_to_SSE_16x8i, rec(
    isa_from := SKLR_64x1i,
    isa_to   := SSE_16x8i(T_Int(8)),

    code     := (self, y, x, opts) >> let(
                    xx := (offs) -> vtref(self.isa_from.t, x, offs),
                    yy := (offs) -> vtref(self.isa_to.t,   y, offs),
                    a  := var.fresh_t("U", T_UInt(64)),
                    b0 := var.fresh_t("U", T_UInt(64)),
                    b1 := var.fresh_t("U", T_UInt(64)),
                    mask := T_UInt(64).value(1 + 2^8 + 2^16 + 2^24 + 2^32 + 2^40 + 2^48 + 2^56),
                    decl([a,b0,b1], chain(
                        assign( a, tcast(a.t, xx(0)) ),
                        assign( b0, bin_and(            a, mask)),
                        assign( b1, bin_and(bin_shr(a, 1), mask)),
                        assign( yy(0), tcast(self.isa_to.t, vpack(b0, b1))),
                        assign( b0, bin_and(bin_shr(a, 2), mask)),
                        assign( b1, bin_and(bin_shr(a, 3), mask)),
                        assign( yy(1), tcast(self.isa_to.t, vpack(b0, b1))),
                        assign( b0, bin_and(bin_shr(a, 4), mask)),
                        assign( b1, bin_and(bin_shr(a, 5), mask)),
                        assign( yy(2), tcast(self.isa_to.t, vpack(b0, b1))),
                        assign( b0, bin_and(bin_shr(a, 6), mask)),
                        assign( b1, bin_and(bin_shr(a, 7), mask)),
                        assign( yy(3), tcast(self.isa_to.t, vpack(b0, b1)))
                    )))
));

Class(SKLR_32x1i_to_SSE_4x32f_f32, SKLR_32x1i_to_SSE_16x8i, rec(
    isa_from := SKLR_32x1i,
    isa_to   := SSE_4x32f(T_Real(32)),

    code     := (self, y, x, opts) >> let(
                    xx := (offs) -> vtref(self.isa_from.t, x, offs),
                    yy := (offs) -> vtref(self.isa_to.t,   y, offs),
                    ti := TVect(T_Int(32), 4),
                    tf := TVect(T_Real(32), 4),
                    a  := var.fresh_t("U", T_UInt(32)),
                    b  := var.fresh_t("U", ti),
                    decl( [a, b], chain(
                        assign( a, tcast(a.t, xx(0)) ),
                        assign( b, vpack(a, bin_shr(a, 8), bin_shr(a, 16), bin_shr(a, 24)) ),
                        assign(yy(0), tcast(tf, bin_and( tf.one(), arith_shr(bin_shl(b, 31), 31) ))),
                        assign(yy(1), tcast(tf, bin_and( tf.one(), arith_shr(bin_shl(b, 30), 31) ))),
                        assign(yy(2), tcast(tf, bin_and( tf.one(), arith_shr(bin_shl(b, 29), 31) ))),
                        assign(yy(3), tcast(tf, bin_and( tf.one(), arith_shr(bin_shl(b, 28), 31) ))),
                        assign(yy(4), tcast(tf, bin_and( tf.one(), arith_shr(bin_shl(b, 27), 31) ))),
                        assign(yy(5), tcast(tf, bin_and( tf.one(), arith_shr(bin_shl(b, 26), 31) ))),
                        assign(yy(6), tcast(tf, bin_and( tf.one(), arith_shr(bin_shl(b, 25), 31) ))),
                        assign(yy(7), tcast(tf, bin_and( tf.one(), arith_shr(bin_shl(b, 24), 31) )))
                    ))
                )
));

Class(SKLR_64x1i_to_SSE_4x32f_f32, SKLR_32x1i_to_SSE_16x8i, rec(
    isa_from := SKLR_64x1i,
    isa_to   := SSE_4x32f(T_Real(32)),

    code     := (self, y, x, opts) >> let(
                    xx := (offs) -> vtref(self.isa_from.t, x, offs),
                    yy := (offs) -> vtref(self.isa_to.t,   y, offs),
                    ti := TVect(T_Int(32), 4),
                    tf := TVect(T_Real(32), 4),
                    a  := var.fresh_t("U", T_UInt(64)),
                    b0 := var.fresh_t("U", ti),
                    b1 := var.fresh_t("U", ti),
                    shift := (t, n) -> tcast(T_Int(32), bin_shr(t, n)),
                    decl( [a, b0, b1], chain(
                        assign( a, tcast(a.t, xx(0)) ),
                        assign( b0, vpack(shift(a,  0), shift(a,  8), shift(a, 16), shift(a, 24)) ),
                        assign( b1, vpack(shift(a, 32), shift(a, 40), shift(a, 48), shift(a, 56)) ),
                        assign(yy( 0), tcast(tf, bin_and( tf.one(), arith_shr(bin_shl(b0, 31), 31) ))),
                        assign(yy( 1), tcast(tf, bin_and( tf.one(), arith_shr(bin_shl(b1, 31), 31) ))),
                        assign(yy( 2), tcast(tf, bin_and( tf.one(), arith_shr(bin_shl(b0, 30), 31) ))),
                        assign(yy( 3), tcast(tf, bin_and( tf.one(), arith_shr(bin_shl(b1, 30), 31) ))),
                        assign(yy( 4), tcast(tf, bin_and( tf.one(), arith_shr(bin_shl(b0, 29), 31) ))),
                        assign(yy( 5), tcast(tf, bin_and( tf.one(), arith_shr(bin_shl(b1, 29), 31) ))),
                        assign(yy( 6), tcast(tf, bin_and( tf.one(), arith_shr(bin_shl(b0, 28), 31) ))),
                        assign(yy( 7), tcast(tf, bin_and( tf.one(), arith_shr(bin_shl(b1, 28), 31) ))),
                        assign(yy( 8), tcast(tf, bin_and( tf.one(), arith_shr(bin_shl(b0, 27), 31) ))),
                        assign(yy( 9), tcast(tf, bin_and( tf.one(), arith_shr(bin_shl(b1, 27), 31) ))),
                        assign(yy(10), tcast(tf, bin_and( tf.one(), arith_shr(bin_shl(b0, 26), 31) ))),
                        assign(yy(11), tcast(tf, bin_and( tf.one(), arith_shr(bin_shl(b1, 26), 31) ))),
                        assign(yy(12), tcast(tf, bin_and( tf.one(), arith_shr(bin_shl(b0, 25), 31) ))),
                        assign(yy(13), tcast(tf, bin_and( tf.one(), arith_shr(bin_shl(b1, 25), 31) ))),
                        assign(yy(14), tcast(tf, bin_and( tf.one(), arith_shr(bin_shl(b0, 24), 31) ))),
                        assign(yy(15), tcast(tf, bin_and( tf.one(), arith_shr(bin_shl(b1, 24), 31) )))
                    ))
                )
));

# NOTE: assumption that most significant bit is set for non zero numbers:
#
# SSE_16x8i_i8_to_SKLR_32x1i can be implemented as
#    neg(vmovemask_16x8i(eq(self.isa_from.t.zero(), xx(2*i))))
# and later simplified if xx(2*i) comes from comparision, yet we cannot match this situation
# because it's unlikely that we will have comparision propagated into this expression.
# Another way is to make T_Bool and simplify expression above by looking at data type.

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

ISA_Bridge.add(Class(CVT_SKLR_16x1i_SSE_16x8ui, CVT_SKLR_16x1i_SSE_16x8i, rec(
    isa_from := SSE_16x8i(T_UInt(8)),
)));

Class(SKLR_32f_to_SKLR_32x1i, ISA_Bridge_I, rec(
    isa_from    := SKLR(T_Real(32)),
    isa_to      := SKLR_32x1i,
    granularity := self >> self.isa_to.v,

    code := (self, y, x, opts) >> let(
                j  := Ind(self.isa_to.v),
                xt := self.isa_from.t,
                yt := T_UInt(self.isa_to.v),
                yy := (offs) -> vtref(self.isa_to.t, y, offs),
                a  := var.fresh_t("U", yt),
                decl([a], chain(
                    assign(a, a.t.zero()),
                    loop(j, j.range,
                        assign(a, bin_or(a, cond(eq(nth(x, j), xt.zero()), yt.zero(), bin_shl(yt.one(), j))))),
                    assign(yy(0), a)
                ))
            ),
));

Class(SKLR_32f_to_SKLR_64x1i, SKLR_32f_to_SKLR_32x1i, rec(
    isa_from := SKLR(T_Real(32)),
    isa_to   := SKLR_64x1i
));

Class(SKLR_64x1i_to_SKLR_32f, ISA_Bridge_I, rec(
    isa_from    := SKLR_64x1i,
    isa_to      := SKLR(T_Real(32)),
    granularity := self >> self.isa_from.v,

    code := (self, y, x, opts) >> let(
                j  := Ind(self.isa_from.v),
                xx := (offs) -> vtref(self.isa_from.t, x, offs),
                ti := T_UInt(self.isa_from.v),
                tf := self.isa_to.t,
                a  := var.fresh_t("U", ti),
                decl( [a], chain(
                    assign( a, tcast(a.t, xx(0)) ),
                    loop(j, j.range,
                        assign( nth(y, j), tcvt( tf, bin_and(bin_shr(a, j), a.t.one())))
                    ).unroll()
                ))
            ),
));

Class(SKLR_32x1i_to_SKLR_32f, SKLR_64x1i_to_SKLR_32f, rec(
    isa_from := SKLR_32x1i,
    isa_to   := SKLR(T_Real(32))
));

Class(SKLR_64x1i_to_SKLR_8i, SKLR_64x1i_to_SKLR_32f, rec(
    isa_from := SKLR_64x1i,
    isa_to   := SKLR(T_Int(8))
));

Class(SKLR_32x1i_to_SKLR_8i, SKLR_64x1i_to_SKLR_32f, rec(
    isa_from := SKLR_32x1i,
    isa_to   := SKLR(T_Int(8))
));

ISA_Bridge.add(Class(CVT_SKLR_32x1i_SKLR_16x1i, ISA_Bridge_I, rec(
    isa_from    := SKLR_16x1i,
    isa_to      := SKLR_32x1i,
    code := (self, y, x, opts) >>
        assign(self._y(y,0), bin_or(tcvt(T_UInt(32), self._x(x,0)), bin_shl(tcvt(T_UInt(32), self._x(x,1)), 16)) ),
)));

ISA_Bridge.add(Class(CVT_SKLR_64x1i_SKLR_16x1i, ISA_Bridge_I, rec(
    isa_from    := SKLR_16x1i,
    isa_to      := SKLR_64x1i,
    code := (self, y, x, opts) >>
        assign(self._y(y,0), bin_or(
                     tcvt(T_UInt(64), self._x(x,0)),
             bin_shl(tcvt(T_UInt(64), self._x(x,1)), 16),
             bin_shl(tcvt(T_UInt(64), self._x(x,2)), 32),
             bin_shl(tcvt(T_UInt(64), self._x(x,3)), 48)
        )),
)));


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

# Non-Terminals (Transforms)
# ==========================
# MP, from 03/24/01, GAPv3.4.4

#F Non-Terminaltable (Transform Table)
#F -----------------------------------
#F

#F The non-terminal table represents all symbols that can be used
#F to define an spl of type "nonTerminal". The table is meant to be
#F extended. In contrast to spls of type "symbol", an spl of type
#F "nonTerminal" has no equivalent in the SPL language. The non-terminals
#F usually represent discrete signal transforms and are used by the
#F formula generator to represent non-terminal expressions in SPL objects.
#F These are meant to be further expanded. Thus, an SPL object can be exported
#F to a valid SPL program iff it contains no non-terminals.
#F Every non-terminal is represented by a record contained in the
#F following list NonTerminalTableSPL.
#F The record has the form:
#F rec(
#F   isNonTerminal,   true, identifying a non-terminal
#F   NonTerminalOps,  operations record
#F   symbol,          a string for the symbol
#F   CheckParams,     a function for checking and canonifying
#F                    of the parameters
#F   Dim,             a function to return the dimension from
#F                    the parameters in canonical form
#F   Params,          a function that returns a shorter representation
#F                    of parameters used for printing in gap
#F   Terminate SPL,   a function to convert into an spl without non-terminals
#F   Transposed,      a function for transposing
#F   isRealTransform, = true if the transform is considered as real transform,
#F                    i.e., the matrix is real for all choices of parameters
#F )
#F
#F the following fields are optional:
#F for verification:
#F SmallRandom        a function that returns a parameter choice for a small
#F                    instantiation of the transform, i.e., one that can be verified
#F                    against the definition
#F LargeRandom        correspondingly, a function that returns a parameter choice
#F                    for a large instantiation of the transform, used for verifying
#F                    on a random vector
#F
#F HashIndexSPL       a function that returns index used for hashing. It is stored
#F                    in .hashIndex field. Overrides the standard
#F                    HashIndex function used in the search.
#F
#F hashIndex          overrides HashIndex function used in
#F                    search. Precomputed by optional HashIndexSPL
#F                    function and should not be specified manually.
#F
#F
#F Nonterminals are created in the directory formgen/transforms
#F
#F The Non-Terminal table can easily be extended by adding new non-terminals
#F (transforms).
#F

#F Non-Terminal Table
#F ------------------
#F

#F NonTerminalTableSPL
#F   is the set of all known non-terminals (transforms)
#F
NonTerminalTableSPL := Set ( [ ] );

#F NonTerminalListSPL
#F   is a set containing all known non-terminals.
#F   Note that NonTerminalListSPl is ordered exactly as NonTerminalTableSPL.
#F
NonTerminalListSPL := Set( [ ] );

#F Adding new Non-Terminals (Transforms)
#F -------------------------------------
#F

#F AddNonTerminal ( <non-terminal> )
#F   adds <non-terminal> to the non-terminal table. This function is
#F   meant to import new non-terminals into NonTerminalTableSPL,
#F   and NonTerminalListSPL, and update .index field.
#F
#F Note: Currently, this function is invoked automatically, when one
#F       for example, defines new breakdown rules for a non-terminal,
#F
AddNonTerminal := function ( S )
  local i;
  Constraint(IsNonTerminal(S));
  if not IsBound(S.__bases__) then
      Error("Old-style non-terminals are not supported anymore");
  fi;

  if S in NonTerminalTableSPL then RemoveSet(NonTerminalTableSPL, S); fi;
  AddSet(NonTerminalTableSPL, S);
  AddSet(NonTerminalListSPL, S.name);
  for i in [1..Length(NonTerminalTableSPL)] do
      NonTerminalTableSPL[i].index := i;
  od;
  #UpdateApplicableTable();
end;


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

Declare(StreamGath);
Declare(StreamScat);
Declare(STensor);
Declare(SIterDirectSum);

Import(compiler);

Class(RootsOfUnity, Sym, rec(
    abbrevs := [ (n) -> [n]],

    def := (n) -> let(j := Ind(n),
	Diag(Lambda(j, omega(n, j)).setRange(Complexes))),

    tolist := self >> self.lambda().tolist(),
    lambda := self >> self.obj.element,
    sums := self >> Diag(self),
    isFunction := true,
    range := self >> TComplex,
    domain := self >> self.params[1],

    transpose := self >> self,

    isReal := self >> (self.params[1] <= 2),

    exportParams := self >> self.params

));

Class(streamDiagFunc, FuncClass, rec(
    def := (N, n, r) -> rec(
        N := N/2,
        n := N,
	r := r,
    ),

    tolist := self >> self.lambda().tolist(),
    lambda := self >> let(
	i := Ind(self.n),
	Lambda(i, imod(i, self.params[2]) * idiv(i, self.params[2]))),
   domain := self >> self.n,
   range := self >> self.N
));

# ==========================================================================
# StreamGath(<i>, <range>, <bs>) - vector gather (read) matrix
# ==========================================================================
Class(StreamGath, BaseMat, SumsBase, rec(
    #-----------------------------------------------------------------------
#    rChildren := self >> [self.func],
#    rSetChild := rSetChildFields("func"),
    rChildren := self >> [self.i],
    rSetChild := rSetChildFields("i"),
    #-----------------------------------------------------------------------
    new := (self, i, range, bs) >> SPL(WithBases(self,
        rec(dimensions := [bs, bs*range], bs := bs, i := i, range := range))),
    #-----------------------------------------------------------------------
    sums := self >> self,
    area := self >> self.bs * self.range,
    isReal := self >> true,
    transpose := self >> StreamScat(self.i, self.range, self.bs),
    #-----------------------------------------------------------------------
    print := (self,i,is) >> Print(self.name, "(", self.i, ", ", self.range,", ", self.bs,")"),
    #-----------------------------------------------------------------------
    toAMat := self >> Gath(fTensor(fBase(self.range, self.i), fId(self.bs))).toAMat()
    # make verilog
));

# ==========================================================================
# StreamScat(<func>, <v>) - vector scatter (write) matrix
# ==========================================================================
Class(StreamScat, BaseMat, SumsBase, rec(
    #-----------------------------------------------------------------------
    rChildren := self >> [self.i],
    rSetChild := rSetChildFields("i"),
    #-----------------------------------------------------------------------
    new := (self, i, range, bs) >> SPL(WithBases(self,
        rec(dimensions := [bs*range, bs], bs := bs, i := i, range := range))),
    #-----------------------------------------------------------------------
    sums := self >> self,
    area := self >> self.func.domain() * self.v,
    isReal := self >> true,
    transpose := self >> StreamGath(self.i, self.range, self.bs),
    #-----------------------------------------------------------------------
    print := (self,i,is) >> Print(self.name, "(", self.i, ", ", self.range,", ", self.bs,")"),
    #-----------------------------------------------------------------------
    toAMat := self >> TransposedAMat(StreamGath(self.i, self.range, self.bs).toAMat())
#    code := (self, y, x) >> VTensor(Scat(self.func), self.v).code(y,x)
));

RCExpandData := function(value)
   return ConcatList(value.v, x->[ReComplex(Complex(x.v)), ImComplex(Complex(x.v))]);
end;

#FixLocalize := (o) -> SubstBottomUpRules(o, [
#    [[Diag, @(1)], e -> Diag(GenerateData(ResolveMemos(e.element))), "FixLocalize" ]]);

Declare(BRAMPerm);

# ===============================================
# CodeBlock(<SPL>)
# ==============================================
Class(CodeBlock, BaseContainer, SumsBase, rec(
    isBlock:=true,

    sums := self >> CopyFields(self, rec(_children := [self.child(1).sums()])),

    code := false,

	createCode := meth(self)
	   local i, s, res, d, dl;
	   if (IsBound(self._children[1].isBlock)) then
	     res := self._children[1].createCode();
	   else
	     s := self.child(1);
	     res := CopyFields(self, rec(code := CodeSums(s, StreamDefaults)));
	   fi;
       res.code := BinSplit(res.code);
	   return res;
	end,

    rChildren := self >> [self._children[1], self.code],

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

    from_rChildren := (self, rch) >> CopyFields(ObjId(self)(rch[1]), rec(code:=rch[2])),

    print := (self,i,is) >> Print(self.name, "(", Cond(self.code<>false, self.code.print(i,is), self._children[1]), ")")
));

Class(StreamSum, ISum, BaseIterative, SumsBase);

# A streaming I tensor J primitive
Class(StreamIxJ, BaseMat, SumsBase, rec(
    new := (self, k, p, bs) >> SPL(WithBases(self, rec(dimensions:=[k*p,k*p], element:=k, p:=p, bs:=bs))),
    print := (self,i,is) >> Print(self.name, "(", self.element, ", ", self.p, ", ", self.bs, ")"),
    toAMat := self >> let(J := When(self.element = 2, I(2), Compose(Tensor(I(2), L(self.element/2, self.element/4)), L(self.element, 2))), Tensor(I(self.p), J).toAMat()),
    transpose := self >> self,
    area := self >> self.element,
    sums := self >> self,
    rChildren := self >> [],
    dimensions := self >> [self.element * self.p, self.element * self.p],
    dims := self >> self.dimensions,
   isReal := True
));

# parameters:
#   L[1]: SPL
#   L[2]: p: p in I_p x SPL
#   L[3]: bs: stream buffer size
Class(STensor, Tensor,
  rec(
    new := (self, L) >> SPL(WithBases(self, rec(
            bs := L[3],
            dimensions := L[1].dimensions * Cond(IsValue(L[2]), L[2].ev(), L[2]),
            p := Cond(L[1].dims()[1] >= L[3], L[2], L[2]*L[1].dims()[1]/L[3]),
            _children := [Checked(L[1].dims()[1] = L[1].dims()[2],
                          Cond(L[1].dims()[1] >= L[3], L[1], Tensor(I(L[3] / L[1].dims()[1]), L[1]))), Cond(L[1].dims()[1] >= L[3], L[2], L[2]*L[1].dims()[1]/L[3]), L[3]]
#            _children := [Checked(L[1].dims()[1] = L[1].dims()[2],
#                          Cond(L[1].dims()[1] >= L[3], L[1], Tensor(I(L[3] / L[1].dims()[1]), L[1])))]
                          ))),

    print := (self,i,is) >> Print(self.name, "(", self.child(1), ", ", self.p, ", ", self.bs, ")"),
    createCode := self >> Cond(IsBound(self.child(1).createCode), STensor(self.child(1).createCode(), self.p, self.bs), self),
    isPermutation := False,
    toAMat := self >> Tensor(I(self.p), self.child(1)).toAMat(),
    sums := self >> self,
    dims := self >> Cond(IsValue(self.p), self.p.ev(), self.p) * self.child(1).dims()
));

# parameters: I_m x L^n_r x I_p
#   L[1]: n
#   L[2]: r
#   L[3]: p
#   L[4]: m
#   L[5]: stream buffer size
Class(SIxLxI, Tensor,
  rec(
    new := (self, l) >> SPL(WithBases(self, rec(
            _children := [],
            dimensions := [l[1] * l[3] * l[4], l[1] * l[3] * l[4]],
            bs := l[5],
            p := l[3],
            m := l[4],

    createCode := self >> self,
    print := (self,i,is) >> Print(self.name, "(", l[1], ", ", l[2], ", ", self.p, ", ", self.m, ", ", self.bs, ")"),
    isPermutation := False,
    toAMat := self >> Tensor(Tensor(I(self.m), L(l[1], l[2]), I(self.p))).toAMat(),
    transpose := self >> self,
    dims := self >> self.dimensions,
#    dimensions := l[1] * l[3] * l[4],
    sums := self >> self,
)))));

# ==========================================================================
# Streaming iterative direct sum
# SIterDirectSum(<var>, <range>, <spl>, <stream buffer size>)
#    Size of the <spl> must not depend on <var>.
# ==========================================================================
Class(SIterDirectSum, BaseIterative, rec(
    abbrevs := [ (v,d,e, bs) -> [v,d,e,bs],
	         (v,e,bs) -> [v, v.range, e, bs]],
    new := meth(self, var, domain, expr, bs)
        local res;
        var.isLoopIndex := true;
        res := SPL(WithBases(self, rec(
#		_children := [Checked(expr.dims()[1] <= bs and
		_children := [Checked(expr.dims()[1] = expr.dims()[2],
				      Cond(expr.dims()[1] >= bs,
					  expr,
					  IterDirectSum(Ind(bs/expr.dims()[1]), expr)))],
		var := var,
		bs := bs,
		dimensions := expr.dimensions * domain,
		domain := domain)));
        return res;
        end,

    rChildren := self >> self._children,
    rSetChild := meth(self, n, what)  self._children[n] := what; end,
    child := (self, n) >> self._children[n],
    children := self >> self._children,

#         Constraint(IsSPL(expr));
# 	Constraint(not IsList(domain));
# 	#domain := toRange(domain);
# 	if domain = 1 then
# 	    return SubstBottomUp(expr, @.cond(e->Same(e,var)), e->V(0));
# 	fi;
# 	var.range := domain;
#         return SPL(WithBases( self,
# 	    rec( expr := expr,
# 		 var := var,
# 		 domain := domain,
# 		 dimensions := domain * Dimensions(expr) )));
#     end,

    print := (self, i, is) >> Print(self.name, "(", self.var, ", ",  self.domain, ", ", self.child(1).print(i+is, is), ", ", self.bs, ")"),

    toAMat := self >> self.unroll().toAMat(),

    unroll := self >> DirectSum(
	List([ 0 .. self.domain - 1 ],
	    index_value -> SubstBottomUp(
		Copy(self._children[1]), self.var,
		e -> V(index_value)))),

    dims := self >> self.dimensions,
    #-----------------------------------------------------------------------
    isPermutation := False,
    #-----------------------------------------------------------------------
    createCode := self >> Cond(IsBound(self._children[1].createCode),
	SIterDirectSum(self.var, self.domain, self._children[1].createCode(), self.bs), self),

    transpose :=   # we use inherit to copy all fields of self
        self >> Inherit(self, rec(
		   expr := TransposedSPL(self._children[1]),
		   dimensions := Reversed(self.dimensions))),

    sums := self >> self
));

NoPull.createCode := meth ( self )
    local i, s, res, d, dl;
    if IsBound(self._children[1].isBlock)  then
        res := self._children[1].createCode();
    else
        s := self.child(1);
        res := SPL(CopyFields(self, rec(
                code := CodeSums(self.dims()[1] * 2, s) )));
    fi;
    for i  in [ 1 .. 10 ]  do
        BinSplit(res.code);
    od;
    return res;
end;

Class(Omega, BaseIterative, rec(

    abbrevs := [ (N,n,i) ->[N, n, i] ],

    new := (self, N, n, i) >> let(j := Ind(n),
        Lambda(j, omega(N, i*j)).setRange(TComplex)),

    rChildren := self >> self._children,
    rSetChild := meth(self, n, what)  self._children[n] := what; end,
    child := (self, n) >> self._children[1],
    children := self >> self._children,

    tolist := self >> self.lambda().tolist(),
    lambda := self >> self.obj.element,
    sums := self >> Diag(self),
    isFunction := true,
    range := self >> TComplex,
    domain := self >> self.params[1],

    transpose := self >> self,

#    isReal := self >> (self.params[1] = self.params[2]  or self.params[2] = 1),

    exportParams := self >> self.params
));

Prm.createCode := (self) >> SPL(CopyFields(self, rec(code := CodeSums(self, StreamDefaults))));
ISum.createCode := (self) >> CodeBlock(self).createCode();

Tensor.createCode := self >>
   Tensor(
      List([1..Length(self.children())], i->
        Cond(IsBound(self.child(i).createCode),
	    self.child(i).createCode(),
	    self.child(i)
         )
      )
   );

Tw1Stride := (n,d,k,which,stride) -> Checked(IsPosIntSym(n), IsPosIntSym(d), IsPosIntSym(k),
    fCompose(dOmega(n,k), diagTensor(dLin(n/d, 1, 0, TInt), dLin(d/stride, 1, which, TInt))));

# fReal := (f) -> fCompose(RCData(f), fTensor(fId(f.n), fBase(2, 0)));
Class(fReal, FuncClass, rec(
   def := (f) -> rec(N := f.range(), n := f.domain(), f := f),
   tolist := self >> self.lambda().tolist(),
   lambda := self >> fCompose(RCData(self.f), fTensor(fId(self.n), fBase(2,0))).lambda(),
   domain := self >> self.n,
   range := self >> self.N
));

# fImag := (f) -> fCompose(RCData(f), fTensor(fId(f.n), fBase(2, 1)));
Class(fImag, FuncClass, rec(
   def := (f) -> rec(N := f.range(), n := f.domain(), f := f),
   tolist := self >> self.lambda().tolist(),
   lambda := self >> fCompose(RCData(self.f), fTensor(fId(self.n), fBase(2,1))).lambda(),
   domain := self >> self.n,
   range := self >> self.N
));

Class(Stream1DFT2, BaseMat, SumsBase, rec(
    new := (self) >> SPL(WithBases(self, rec(dimensions:=[2,2]))),
    print := (self,i,is) >> Print(self.name),
    toAMat := self >> DFT(2,1).toAMat(),
    transpose := self >> self,
    sums := self >> self,
    rChildren := self >> [],
    dimensions := [2,2],
    dims := self >> self.dimensions,
));

Class(StreamPadDiscard, BaseMat, SumsBase, rec(
    abbrevs := [(n, m, w_in, w_out) -> [n, m, w_in, w_out]],

    new := (self, n, m, w_in, w_out) >> SPL(WithBases(self, rec(
          _children := [n,m, w_in, w_out],
          dimensions := [n, m],
    ))),

    rChildren := self >> self._children,
	rSetChild := meth(self, n, what) self._children[n] := what; end,
	child     := (self, n) >> self._children[n],
	children  := self >> self._children,

    toAMat := self >> RI(self.child(1), self.child(2)).toAMat(),
    dims := self >> [self.child(1), self.child(2)],

));


text:
# Several ways to do it
"Goodbye, World!";

Print("Goodbye, World!\n"); # No EOL appended

Display("Goodbye, World!");

f := OutputTextUser();
WriteLine(f, "Goodbye, World!\n");
CloseStream(f);

max_stars_repo_name: djgoku/RosettaCodeData

max_issues_repo_name: djgoku/RosettaCodeData

max_forks_repo_name: djgoku/RosettaCodeData


text:
IsPrimeTrial := function(n)
   local k, m;
   if n < 5 then
      return (n = 2) or (n = 3);
   fi;
   if RemInt(n, 2) = 0 then
      return false;
   fi;
   m := RootInt(n);
   k := 3;
   while k <= m do
      if RemInt(n, k) = 0 then
         return false;
      fi;
      k := k + 2;
   od;
   return true;
end;

Filtered([1 .. 100], IsPrimeTrial);
# [ 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97 ]


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

#F Circulant(<n>, <filt-func>, <valuation>) -- n x n  circulant matrix non-terminal
#F
#F The circulant matrix with coefficients of the first row given by
#F the filter function <filt-func>.
#F
#F <valuation> is the offset, default = 0, meaning no offset.
#F
#F Examples:
#F
#F spiral> PrintMat(MatSPL( Circulant(5, FList(TInt, [1..3]), 0)));
#F [ [ 1, 2, 3,  ,   ],
#F   [  , 1, 2, 3,   ],
#F   [  ,  , 1, 2, 3 ],
#F   [ 3,  ,  , 1, 2 ],
#F   [ 2, 3,  ,  , 1 ] ]
#F spiral> PrintMat(MatSPL( Circulant(5, FList(TInt, [1..3]), -1)));
#F [ [ 2, 3,  ,  , 1 ],
#F   [ 1, 2, 3,  ,   ],
#F   [  , 1, 2, 3,   ],
#F   [  ,  , 1, 2, 3 ],
#F   [ 3,  ,  , 1, 2 ] ]

Class(Circulant, NonTerminal, rec(
    _short_print := true,

    abbrevs := [
     (L)     -> [ toSize(L), toFunc(L), toValuation(L) ],
     (n,L)   -> [ Checked(IsPosInt(n), n), toFunc(L), toValuation(L) ],
     (n,L,v) -> [ Checked(IsPosInt(n), n), toFunc(L), Checked(IsInt(v), v) ],
    ],

    dims := self >> Replicate(2, self.params[1]),

    setData := meth(self, newdata) self.params[2] := newdata; return self; end,

    column := self >> let(
	coeffs := List(self.params[2].tolist(), x->x.ev()),
	valuation := self.params[3],
	CircularWrap([coeffs, valuation], self.params[1])),

    terminate := meth(self)
      local i, j, n, L, Ls, mat;
      L := self.column();
      n := Length(L);
      mat := [ ];
      for i in [0..n-1] do
          Ls := Replicate(n, 0);
	  for j in [0..n-1] do Ls[(i-j) mod n + 1] := L[j+1]; od;
	  Add(mat,Ls);
      od;
      return Mat(mat);
    end,

    filtlen := self >> self.params[2].domain(),

    transpose := self >> Circulant(
	self.params[1],
	fCompose(self.params[2], J(self.filtlen())),
	-self.params[3]-self.filtlen()+1),

    isReal := self >> not IsComplexT(self.params[2].range()),

    hashAs := self >> let(
	t := ObjId(self)(
	    self.params[1], fUnk(self.params[2].range(), self.params[2].domain()), self.params[3]),
	When(self.transposed, t.transpose(), t)),

    SmallRandom := () -> let(n := Random([2..16]),
	List([1..n], i -> Float(Random([-2^10..2^10]), Random([-20..10])))),

    LargeRandom := () -> let(n := Random([16, 32, 64]),
	List([1..n], i -> Float(Random([-2^10..2^10]), Random([-20..10]))))
));

##
## Rules
##
RulesFor(Circulant, rec(
    ###################################################################
    ## Time domain methods
    ###################################################################

    #F Circulant_Base: (base case)
    #F
    Circulant_Base := rec (
	info             := "Circulant -> Mat",
	forTransposition := false,
	isApplicable     := P -> P[1] <= 32,
	allChildren      := P -> [[ ]],
	rule := (P, C) -> ApplyFunc(Circulant, P).terminate()
    ),

#     Circulant_toFilt := rec (
# 	info             := "Circulant -> VStack(.., Filt, ...)",
# 	forTransposition := false,
# 	isApplicable     := P -> P[1] > 2 and (P[1]-P[2].domain()) > 2 and P[2].domain() <= 32 and P[3] <= 0,
# 	allChildren      := P -> let(
# 	    n := P[1], l := -P[3], r := P[3] + P[2].domain() - 1,
# 	    [[ Filt(n-l-r, Poly(List(P[2].tolist(),EvalScalar), P[3])) ]]),

# 	rule := (P, C) -> let(
# 	    coeffs := P[2], nc := coeffs.domain(), n := P[1], l := -P[3],
# 	    r := -l + nc - 1,
# 	    j := Ind(l), k := Ind(r), v := Ind(nc),
# 	    VStack(
# 		BB(ISum(j, j.range,
# 		    Scat(fBase(j)) *
# 		    RowVec(fCompose(coeffs, Lambda(v, cond(leq(v,r+j), v-j+l, v-j+l-nc))))) *
# 		    Gath(n, nc,	Lambda(v, cond(leq(v,r+j), v, v+(n-l-r)-1 )))),
# 		C[1],
# 		BB(ISum(k, k.range,
# 		    Scat(fBase(k)) *
# 		    RowVec(fCompose(coeffs, Lambda(v, cond(leq(v,k), v-k+(nc-1), v-k-1)))) *
# 		    Gath(n, nc, Lambda(v, cond(leq(v,k), v, v+(n-l-r)-1 )))))))
#     ),

    #F Circulant_Blocking
    #F
    #F Circulant -> Blocks (Toeplitzes)
    #F
    #F Circulant is partitioned into square blocks of sizes
    #F given by all proper divisors of the input size.
    #F
    #F NOTE: this rule is invalid, Circulant_BlockingDense seems valid
    #F        whats the difference between the two ???
    Circulant_Blocking := rec(
	info             := "Circulant -> Toeplitz",
	forTransposition := false,
	isApplicable     := P -> let(n := P[1],
	    n > 2 and not IsPrime(n)), # and n <= 2*Length(P[2])),

	allChildren := P -> let(
	    divs := DivisorPairs(P[1]),
	    ratio := P[2].domain() / P[1],
	    List(divs, d -> [ Toeplitz(FUnk(2*d[1]-1), 0, CeilingRat(ratio * (2*d[1]-1))) ])
	),

        prepData := (self, L, b, nc) >> let(
	    n := Length(L),
	    Concatenation(
		List([1..nc], k ->
		    List([1..2*b-1], j -> L[((n-k*b+j) mod n) + 1])))),

	rule := (self,P,C,Nonterms) >> let(
	    n    := P[1],
	    b    := Rows(C[1]),
	    nc   := n/b,
	    new_data := FData(self.prepData(ApplyFunc(Circulant, P).column(), b, nc)),
	    i := Ind(nc),
	    k := Ind(nc),
	    ofs := DataInd(TInt, nc * (2*b-1)).setAttr("live_out"),
	    bnum := DataInd(TInt, nc).setAttr("live_out"),
	    gfunc := fCompose(new_data, fAdd(nc*(2*b-1), 2*b-1, ofs)),

	    len := P[2].domain(),
	    val := P[3] mod n,
	    btoep_len := 2*b-1,

	    left := val,
	    right := When(len+b >= n, val-1, (val + (len-1) + (b-1)) mod n),
	    exact_right := (val + (len-1)) mod n,

	    # we "quantize" the condition to tell us which block #s (ofs) are 0
	    lblock := btoep_len * Int(left/b),
	    rblock := btoep_len * Int(right/b),

	    rproj := b*bnum + (b - 1 - left),
	    lproj := b*bnum + (b - 1 - exact_right),
	    nt    := Nonterms[1].setData(gfunc)
	                        .setNonZero(lproj, rproj),

	    open   := leq(rblock+1, lblock),
	    A := leq(lblock, ofs),
	    B := leq(ofs, rblock),
	    isNonZero := logic_or(logic_and(A, B), logic_and(open, logic_or(A, B))),

	    ISum(i, i.range,
		Scat(fTensor(fBase(i), fId(b))) *
		ISumAcc(k, k.range,
		    Data(bnum, imod(nc+k-i, nc),
			Data(ofs, (2*b-1)*bnum,
			    COND(isNonZero, C[1], O(b)))) *
		    Gath(fTensor(fBase(k), fId(b)))))
        )
        # RowDirectSum( nc*b, List([1..nc], i->
        # ColDirectSum( b, List([1..nc], k-> C[((k-i) mod nc) + 1]))))
    ),

    #F Circulant_BlockingDense
    #F
    #F Circulant -> Blocks (Toeplitzes)
    #F
    #F Circulant is partitioned into square blocks of sizes
    #F given by all proper divisors of the input size.
    #F
    Circulant_BlockingDense := rec(
	info             := "Circulant -> Toeplitz",
	forTransposition := false,
	isApplicable     := P -> let(n := P[1],
	    n > 2 and not IsPrime(n) and n <= 2*P[2].domain()),

	allChildren := P -> let(
	    divs := DivisorPairs(P[1]),
	    ratio := P[2].domain() / P[1],
	    List(divs, d -> [ Toeplitz(FUnk(2*d[1]-1), 0, CeilingRat(ratio * (2*d[1]-1))) ])
	),

        prepData := (self, L, b, nc) >> let(
	    n := Length(L),
	    Concatenation(
		List([1..nc], k ->
		    List([1..2*b-1], j -> L[((n-k*b+j) mod n) + 1])))),

	rule := (self,P,C,Nonterms) >> let(
	    n    := P[1],
	    b    := Rows(C[1]),
	    nc   := n/b,
	    new_data := FData(self.prepData(ApplyFunc(Circulant, P).column(), b, nc)),
	    i := Ind(nc),
	    k := Ind(nc),
	    ofs := DataInd(TInt, nc * (2*b-1)).setAttr("live_out"),
	    bnum := DataInd(TInt, nc).setAttr("live_out"),

	    gfunc := fCompose(new_data, fAdd(nc*(2*b-1), 2*b-1, ofs)),

	    len := P[2].domain(),
	    val := P[3] mod n,
	    left := val,
	    exact_right := (val + (len-1)) mod n,

	    rproj := b*bnum + (b - 1 - left),
	    lproj := b*bnum + (b - 1 - exact_right),
	    nt    := Nonterms[1].setData(gfunc)
	                        .setNonZero(lproj, rproj),

	    ISum(i, i.range,
		Scat(fTensor(fBase(i), fId(b))) *
		ISumAcc(k, k.range,
		    Data(bnum, imod(nc+k-i, nc),
			Data(ofs, (2*b-1)*bnum,
			    C[1])) *
		    Gath(fTensor(fBase(k), fId(b)))))
        )
        # RowDirectSum( nc*b, List([1..nc], i->
        # ColDirectSum( b, List([1..nc], k-> C[((k-i) mod nc) + 1]))))
    ),

    #F Circulant_DiagonalizeStep
    #F
    #F Circulant_n -> (F2 tensor I2)
    #F                (Circulant_n/2 dirsum Toeplitz_n/2)
    #F                (F2 tensor I2)
    #F
    #F Computes circulant matrix through DFT.
    #F h is the first column of the circulant
    #F (params of Circulant nonterminal).
    #F
    Circulant_DiagonalizeStep := rec(
	info             := "Circulant_n -> Circulant_n/2 dirsum Toeplitz_n/2",
	forTransposition := false,
	isApplicable     := P -> let(n := P[1],
	    n > 2 and n mod 2 = 0 and n<=16),

	allChildren := P -> let(n := P[1],
	     [[ Circulant(FUnk(n/2)), Toeplitz(FUnk(n-1)) ]]),

	prepData := function(L)
	    local l, L1, L2, circ, toep;
	    l := Length(L);
	    L1 := L{[1     .. l/2]}; # first half
	    L2 := L{[l/2+1 ..  l ]}; # second half

	    circ := 1/2 * (L1 + L2);
	    toep := 1/2 * Concat(Drop(L2,1) - Drop(L1,1),
		                      L1    - L2);
	    return [circ, toep];
	end,

	rule := (self, P, C, Nonterms) >> let(
	    n := P[1],
	    data := self.prepData(ApplyFunc(Circulant,P).column()),
	    circ := Nonterms[1].setData(FData(data[1])),
	    toep := Nonterms[2].setData(FData(data[2])),

	    Tensor(F(2), I(n/2)) *
	    DirectSum(C[1], C[2]) *
	    Tensor(F(2), I(n/2))
	)
    ),

    ###################################################################
    ## Frequency domain methods
    ###################################################################

    #F Circulant_DFT
    #F
    #F Circulant -> invDFT * diag(DFT(h)) * DFT
    #F
    #F Computes circulant matrix through DFT.
    #F h is the first column of the circulant
    #F (params of Circulant nonterminal).
    #F
    Circulant_DFT := rec(
	info             := "Circulant -> invDFT * diag(DFT(h)) * DFT",
	forTransposition := false,
	switch           := false,
	isApplicable     := P -> P[1] > 1,
	allChildren      := P -> [[ transforms.DFT(P[1], -1), transforms.DFT(P[1], 1) ]],

	rule := function(P, C)
            local l, coef, L;
	    L := ApplyFunc(Circulant, P).column();
	    l := Length(L);
	    coef := List([1..l], i -> ComplexAny(L[i]));
	    coef := 1/l * ComplexFFT(coef);
	    return C[1] * Diag(FData(coef)) * C[2];
	end
    ),

    #F Circulant_RDFT
    #F
    #F Circulant -> invRDFT * block(RDFT(h)) * RDFT
    #F
    #F Computes circulant matrix through the Real DFT.
    #F h is the first column of the circulant
    #F (params of nonterminal Circulant).
    #F
    Circulant_RDFT := rec(
	info             := "Circulant -> invRDFT * block(RDFT(h)) * RDFT",
	forTransposition := false,
	switch           := false,

	isApplicable     := P -> let(n := P[1],
	    n > 2 and n <= 64 and (n mod 2 = 0) and
	    ForAll(Factors(n), f -> f in [2, 3])),

	allChildren      := P -> let(n := P[1],
	    [[ RDFT(n).transpose(), RDFT(n) ]]),

	rule := function( P, C )
 	    local i, L, M, X, l, coef, ccoef;

	    L := ApplyFunc(Circulant, P).column();
	    l := Length(L);
	    ccoef := List([1..l], i -> ComplexAny(L[i]));
	    ccoef := ComplexFFT(ccoef);
	    coef  := List([1..Int(l/2)+1], i -> 2/l * ReComplex(ccoef[i]));
	    Append(coef, List([Int(l/2)+2..l], i -> 2/l * (-ImComplex(ccoef[i])) ));

	    X := Diag([ 1/2 * coef[1]]);
	    for i in [2..Int(l/2)] do
	        M := Mat([[coef[i],      coef[l-i+2]],
			  [-coef[l-i+2], coef[i]]]);
		X := DirectSum(X, M);
	    od;
	    if (l mod 2 = 0) then
		X := DirectSum(X, Diag( [1/2 * coef[l/2+1]] ));
	    fi;
	    X := X^perm4(l);
	    return C[1] * X * C[2];
	end
    ),

    #F Circulant_PRDFT
    #F
    #F Circulant -> IPRDFT * diag(PRDFT(h)) * PRDFT
    #F
    #F Computes circulant matrix through the Packed Real DFT.
    #F h is the first column of the circulant
    #F (params of nonterminal Circulant).
    #F
    Circulant_PRDFT := rec(
	info             := " Circulant -> IPRDFT * diag(PRDFT(h)) * PRDFT",
	forTransposition := false,
	switch           := false,

	isApplicable     := P -> true,

	allChildren      := P -> let(n := P[1],
	    [[ IPRDFT(n), PRDFT(n) ]]),

	rule := function( P, C )
 	    local i, n, L, l, coef, D;
	    n := Cols(C[1]);
	    L := ApplyFunc(Circulant, P).column();
	    l := Length(L);
	    coef := List(L, EvalScalar);
	    coef := 1/l * ComplexFFT(coef);
	    D:=Diag(FData(coef{[1..n/2]}));
	    return C[1] * RC(D) * C[2];
	end
    ),

    #F Circulant_DHT
    #F
    #F Circulant -> DHT * block(DHT(h)) * DHT
    #F
    #F Computes circulant matrix through the Discrete Hartley transform.
    #F h is the first column of the circulant
    #F (params of nonterminal Circulant).
    #F
    Circulant_DHT := rec(
	info             := "Circulant -> DHT * block(DHT(h)) * DHT",
	forTransposition := false,
	switch           := false,
	isApplicable     := P -> let(n := P[1], n > 1 and Is2Power(n)),

	allChildren      := P -> let(n := P[1], [[ DHT(n) ]]),

	rule := function(P, C)
	    local i,M,X,l,L,coef,ccoef, rcoef, icoef;
	    L := ApplyFunc(Circulant, P).column();
	    l := Length(L);
	    ccoef := List([1..l], i -> ComplexAny(L[i]));
	    ccoef := ComplexFFT(ccoef);
	    rcoef := List([1..l/2+1], i -> 1/l * ReComplex(ccoef[i]));
	    icoef := List([2..l/2], i -> 1/l * ImComplex(ccoef[i]));
	    X := Diag([rcoef[1]]);
	    for i in [2..Int(l/2)] do
	        M := Mat([[rcoef[i],   -icoef[i-1]],
			  [icoef[i-1], rcoef[i]]]);
		X := DirectSum(X,M);
	    od;
	    X := DirectSum(X, Diag( [rcoef[l/2+1]] ));
	    X := X^perm4(l);
	    return C[1] * X * C[1];
	end
    )
));


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

Class(MultiPtrSumsgenMixin, rec(
    IterHStack :=
        (self, o, opts) >> let(
        	bkcols := Cols(o.child(1)),
        	bkrows := Rows(o.child(1)),
        	nblocks := o.domain,
        	cols := Cols(o), rows := Rows(o),
        	j := o.var,
        	ch := self(o.child(1), opts),

    	    ISum(j, j.range,
        		ch *
        		Gath(fTensor(fBase(nblocks, j), fId(bkcols)))
        	)
        )
));

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 2018-2019, Carnegie Mellon University
# See LICENSE for details

Class(HCOLUnparserMixin, rec(
    copyright := "",
    loopc := (self,o,i,is) >> Checked(IsRange(o.range),
        let(v := o.var, lo := o.range[1], hi := Last(o.range),
            Print(Blanks(i), "for(", v, " = ", lo, "; ", v, self._lt, hi,
            ", ", self(o.cnd, i, is), "; ", v, "++) {\n",
                self(o.cmd,i+is,is),
                Blanks(i), "}\n"))),
    abs := (self, o, i, is) >> self(cond(geq(o.args[1], V(0)), o.args[1], neg(o.args[1])), i, is),
    pow := (self, o, i, is) >> Cond(o.args[1]=2 and o.args[2]>=0, Print("(2 << ",self(o.args[2], i, is), ")"),
     	    	      	     		Print("pow(", self(o.args[1], i, is), ",", self(o.args[2], i, is), ")") ),
    creturn := (self, o, i, is) >> Print(Blanks(i), "return ", self(o.args[1], i, is), ";\n"),
	creturnCond := (self, o, i, is) >> Print(Blanks(i), "if (", self(o.args[1], i, is), ")\n",self(creturn(o.args[2]), i+is, is), "\n"),
    SIMDchain := (self, o, i, is) >> Print(Blanks(i), "#pragma omp simd\n", Blanks(i), "#pragma ivdep\n", self(chain(o.cmds),i,is))
));

Class(HCOLUnparser, HCOLUnparserMixin, CUnparser);
Class(HCOLSSEUnparser, HCOLUnparserMixin, SSEUnparser);

max_stars_repo_name: spiral-software/spiral-package-hcol

max_issues_repo_name: spiral-software/spiral-package-hcol

max_forks_repo_name: spiral-software/spiral-package-hcol