text:
#############################################################################
##
##  ExternalPolymakeCone.gi      JConvex package
##                               Martin Bies
##
##  Copyright 2021               University of Pennsylvania
##
##  A Gap package to do convex geometry by Polymake
##
##  Chapter Cones in Polymake
##
#############################################################################

##############################################################################################
##
##  Section GAP category of PolymakeCones
##
##############################################################################################

DeclareRepresentation( "IsPolymakeConeRep", IsPolymakeCone and IsAttributeStoringRep, [ ] );

BindGlobal( "TheFamilyOfPolymakeCones", NewFamily( "TheFamilyOfPolymakeCones" ) );

BindGlobal( "TheTypeOfPolymakeCone", NewType( TheFamilyOfPolymakeCones, IsPolymakeConeRep ) );

BindGlobal( "JuliaMatrixInt", JuliaEvalString("Matrix{BigInt}") );

BindGlobal( "PolymakeMatrixToGAP", function( mat )
    mat := _Polymake_jl.common.primitive( mat ); # ensure integer coeffs
    return JuliaToGAP( IsList, JuliaMatrixInt( mat ), true );
end );

BindGlobal( "MakePolymakeConeVRep", function( rays, lineality )
    local cone, kwargs;
    cone := Objectify( TheTypeOfPolymakeCone,
            rec( rays := Immutable( rays ),
                 lineality := Immutable( lineality ),
                 number_type := "rational",
                 rep_type := "V-rep" ) );

    kwargs := rec();

    # TODO: verify that kwargs is set correctly, also if one or both lists are empty

    # check for degenerate case
    if Length( rays ) > 0 then
        kwargs.INPUT_RAYS := JuliaMatrixInt( rays );

        # if we also have lineality, add them
        if Length( lineality ) > 0 then
            kwargs.INPUT_LINEALITY := JuliaMatrixInt( lineality );
        fi;
    elif Length( lineality ) > 0 then
        kwargs.INPUT_RAYS := JuliaMatrixInt( lineality );
        kwargs.INPUT_LINEALITY := kwargs.INPUT_RAYS;
    fi;

    cone!.pmobj := CallJuliaFunctionWithKeywordArguments( _Polymake_jl.polytope.Cone, [], kwargs );
    return cone;
end);

BindGlobal( "MakePolymakeConeHRep", function( inequalities, equalities )
    local cone, kwargs;
    cone := Objectify( TheTypeOfPolymakeCone,
            rec( inequalities := Immutable( inequalities ),
                 equalities := Immutable( equalities ),
                 number_type := "rational",
                 rep_type := "H-rep" ) );

    kwargs := rec();

    # TODO: verify that kwargs is set correctly, also if one or both lists are empty

    # check for degenerate case
    if Length( inequalities ) > 0 then
        kwargs.INEQUALITIES := JuliaMatrixInt( inequalities );

        # if we have non-zero equalities, add them
        if Length( equalities ) > 0 then
            kwargs.EQUATIONS := JuliaMatrixInt( equalities );
        fi;
    elif Length( equalities ) > 0 then
        kwargs.INEQUALITIES := JuliaMatrixInt( equalities );
        kwargs.EQUATIONS := kwargs.INEQUALITIES;
    fi;

    cone!.pmobj := CallJuliaFunctionWithKeywordArguments( _Polymake_jl.polytope.Cone, [], kwargs );
    return cone;
end);

##############################################################################################
##
##  Constructors for PolymakeCones
##
##############################################################################################

InstallGlobalFunction( Polymake_ConeByGenerators,
  function( arg )
    local given_rays, i, given_lineality, cone;

    if Length( arg ) = 0 then

        Error( "Wrong input: Please provide some input!" );

    elif Length( arg ) = 1 and IsList( arg[1] ) then

        given_rays := Filtered( arg[ 1 ], row -> not IsZero( row ) );
        if Length( given_rays ) > 0 then

            # non-trivial cone
            return Polymake_ConeByGenerators( arg[ 1 ], [ ] );

        else

            # received the trivial cone
            cone := MakePolymakeConeVRep( [ ], [ ] );
            return cone;

        fi;

    elif Length( arg ) = 2 and IsList( arg[ 1 ] ) and IsList( arg[ 2 ] ) then

        if ( not IsEmpty( arg[ 1 ] ) ) and not ( IsMatrix( arg[ 1 ] ) ) then
            Error( "Wrong input: The first argument should be a Gap matrix!" );
        fi;

        if ( not IsEmpty( arg[ 2 ] ) ) and not ( IsMatrix( arg[ 2 ] ) ) then
            Error( "Wrong input: The second argument should be a Gap matrix!" );
        fi;

        given_rays := Filtered( arg[ 1 ], row -> not IsZero( row ) );
        given_lineality := Filtered( arg[ 2 ], row -> not IsZero( row ) );
        cone := MakePolymakeConeVRep( given_rays, given_lineality );
        return Polymake_CanonicalConeByGenerators( cone );

    fi;

end );

InstallGlobalFunction( Polymake_ConeFromInequalities,
  function( arg )
    local given_ineqs, given_eqs, cone;

    if Length( arg ) = 0 or ForAll( arg, IsEmpty ) then

        Error( "Wrong input: Please provide some input!" );

    elif Length( arg ) = 1 and IsList( arg[ 1 ] ) then

        return Polymake_ConeFromInequalities( arg[ 1 ], [ ] );

    elif Length( arg ) = 2 and IsList( arg[ 1 ] ) and IsList( arg[ 2 ] ) then

        if ( not IsEmpty( arg[ 1 ] ) ) and not ( IsMatrix( arg[ 1 ] ) ) then
            Error( "Wrong input: The first argument should be a Gap matrix!" );
        fi;

        if ( not IsEmpty( arg[ 2 ] ) ) and not ( IsMatrix( arg[ 2 ] ) ) then
            Error( "Wrong input: The second argument should be a Gap matrix!" );
        fi;

        given_ineqs := Filtered( arg[ 1 ], row -> not IsZero( row ) );
        given_eqs := Filtered( arg[ 2 ], row -> not IsZero( row ) );
        cone := MakePolymakeConeHRep( given_ineqs, given_eqs );
        return Polymake_CanonicalConeFromInequalities( cone );

    fi;

end );

##############################################################################################
##
##  Canonicalize cones
##
##############################################################################################

InstallMethod( Polymake_CanonicalConeByGenerators,
               [ IsPolymakeCone ],
  function( cone )
    local rays, lineality;

    if cone!.rep_type = "H-rep" then

        return fail;

    else

        # compute rays
        rays := PolymakeMatrixToGAP( cone!.pmobj.RAYS );

        # extract lineality
        lineality := PolymakeMatrixToGAP( cone!.pmobj.LINEALITY_SPACE );

        # construct the new cone
        return MakePolymakeConeVRep( rays, lineality );

    fi;

end );

InstallMethod( Polymake_CanonicalConeFromInequalities,
               [ IsPolymakeCone ],
  function( cone )
    local ineqs, eqs;

    if cone!.rep_type = "V-rep" then

        return fail;

    else

        # compute facets
        ineqs := PolymakeMatrixToGAP( cone!.pmobj.FACETS );

        # compute linear span
        eqs := PolymakeMatrixToGAP( cone!.pmobj.LINEAR_SPAN );

        # construct the new cone
        return MakePolymakeConeHRep( ineqs, eqs );

    fi;

end );

##############################################################################################
##
##  Conversion of cones
##
##############################################################################################

InstallMethod( Polymake_V_Rep,
               [ IsPolymakeCone ],
  function( cone )
    local rays, lineality;

    if cone!.rep_type = "V-rep" then
        return cone;
    else

        # compute rays
        rays := PolymakeMatrixToGAP( cone!.pmobj.RAYS );

        # extract lineality
        lineality := PolymakeMatrixToGAP( cone!.pmobj.LINEALITY_SPACE );

        # construct the new cone
        return MakePolymakeConeVRep( rays, lineality );

    fi;

end );

InstallMethod( Polymake_H_Rep,
               [ IsPolymakeCone ],
  function( cone )
    local ineqs, eqs;

    if cone!.rep_type = "H-rep" then

        return cone;

    else

        if cone!.rep_type = "V-rep" and cone!.rays = [] then
            return Polymake_ConeFromInequalities( [ [ 0, 1 ], [ -1, -1 ] ] );
        fi;

        # compute facets
        ineqs := PolymakeMatrixToGAP( cone!.pmobj.FACETS );

        # compute linear span
        eqs := PolymakeMatrixToGAP( cone!.pmobj.LINEAR_SPAN );

        # construct the new cone
        return MakePolymakeConeHRep( ineqs, eqs );

    fi;

end );

##############################################################################################
##
##  Attributes of cones
##
##############################################################################################

InstallMethod( Polymake_AmbientSpaceDimension,
              "finding the dimension of the ambient space of the cone",
              [ IsPolymakeCone ],
  function( cone )

    return Length( Polymake_V_Rep( cone )!.rays[1] );

end );

InstallMethod( Polymake_Dimension,
              " returns the dimension of the cone",
            [ IsPolymakeCone ],
  function( cone )

    if Polymake_IsEmpty( cone ) then
        return -1;
    fi;

    return cone!.pmobj.CONE_DIM;

end );

InstallMethod( Polymake_Rays,
              " return the list of generating vertices",
              [ IsPolymakeCone ],
  function( cone )

    return Set( Polymake_V_Rep( cone )!.rays );

end );

InstallMethod( Polymake_Lineality,
              " return the list of generating vertices",
              [ IsPolymakeCone ],
  function( cone )

    return Set( Polymake_V_Rep( cone )!.lineality );

end );

InstallMethod( Polymake_Equalities,
              " return the list of equalities of a cone",
              [ IsPolymakeCone ],
  function( cone )

    return Set( ( Polymake_H_Rep( cone ) )!.equalities );

end );

InstallMethod( Polymake_Inequalities,
              " return the list of inequalities of a cone",
              [ IsPolymakeCone ],
  function( cone )

    return Set( ( Polymake_H_Rep( cone ) )!.inequalities );

end );

InstallMethod( Polymake_RaysInFacets,
              " returns the incident matrix of the rays in the facets",
            [ IsPolymakeCone ],
  function( cone )
    local help_cone, s, res_string, number_rays, ray_list, i, dummy, j, helper;

    help_cone := Polymake_V_Rep( cone );
    number_rays := Length( help_cone!.rays );
    s := JuliaToGAP( IsString, Julia.string( help_cone!.pmobj.RAYS_IN_FACETS ) );

    # process the string
    res_string := SplitString( s, '\n' );
    res_string := List( [ 2 .. Length( res_string ) ], i -> EvalString( res_string[ i ] ) );

    # now construct list of rays -- each facets is a list number_rays entries: 0 indicates that the ray generator is not in the facet and a 1 that it is part of the facet
    ray_list := [];
    for i in [ 1 .. Length( res_string ) ] do
        dummy := List( [ 1 .. number_rays ], i -> 0 );
        for j in [ 1 .. Length( res_string[ i ] ) ] do
            dummy[ res_string[ i ][ j ] ] := 1;
        od;
        Append( ray_list, [ dummy ] );
    od;

    # return the result
    return ray_list;

end );

InstallMethod( Polymake_RaysInFaces,
              " returns the incident matrix of the rays in the faces",
            [ IsPolymakeCone ],
  function( cone )
    local dim, help_cone, rays, rays_in_faces, i, generators, additional_rays, j, converted_additional_rays, pos, pos2, k, help_list;

    # check degenerate case
    dim := Polymake_Dimension( cone );
    if ( dim = 2 ) then

        rays_in_faces := Polymake_RaysInFacets( cone );

    else

        # compute a V-representation of the cone
        help_cone := Polymake_V_Rep( cone );

        # read-off the rays and the list of rays in the facets
        rays := help_cone!.rays;
        rays_in_faces := Polymake_RaysInFacets( help_cone );

        # construct the facets as cones, compute their facets and add them
        # this is thus a recursive function, with (presumably) low performance
        for i in [ 1 .. Length( rays_in_faces ) ] do

            # construct the i-th facet and compute the rays in its facets
            generators := List( Filtered( [ 1 .. Length( rays ) ], j -> not IsZero( rays_in_faces[ i ][ j ] ) ), k -> rays[ k ] );
            additional_rays := RaysInFaces( Cone( generators ) );

            # convert these lists into a list with the rays above
            converted_additional_rays := [];
            for j in [ 1 .. Length( additional_rays ) ] do
                pos := Positions( additional_rays[ j ], 1 );
                pos2 := List( [ 1 .. Length( pos ) ], k -> Position( rays, generators[ pos[ k ] ] ) );
                help_list := List( [ 1 .. Length( rays ) ], k -> 0 );
                for k in pos2 do
                    help_list[ k ] := 1;
                od;
                converted_additional_rays := Concatenation( converted_additional_rays, [ help_list ] );
            od;

            # add this list of rays in the faces of the i-th facet to the above-computed list of rays_in_facets
            rays_in_faces := Concatenation( rays_in_faces, converted_additional_rays  );

        od;

    fi;

    # all all rays, since the overall cone is also a face
    rays_in_faces := Concatenation( rays_in_faces, [ ListWithIdenticalEntries( Length( Polymake_Rays( cone ) ), 1 ) ] );

    # return the result
    return DuplicateFreeList( rays_in_faces );

end );

##############################################################################################
##
##  Properties of PolymakeCones
##
##############################################################################################

InstallMethod( Polymake_IsEmpty,
               "finding if the cone empty is or not",
               [ IsPolymakeCone ],
  function( cone )

    return Length( Polymake_V_Rep( cone )!.rays ) = 0;

end );

InstallMethod( Polymake_IsPointed,
               "finding if the cone is pointed or not",
               [ IsPolymakeCone ],
  function( cone )
    local help_cone;

    help_cone := Polymake_V_Rep( cone );

    return help_cone!.pmobj.POINTED;

end );

##############################################################################################
##
##  Operations with cones
##
##############################################################################################

InstallMethod( Polymake_Intersection,
               "construct intersection of two cones",
               [ IsPolymakeCone, IsPolymakeCone ],
  function( cone1, cone2 )
    local cone1_h, cone2_h, new_ineqs, new_equ;

    cone1_h := Polymake_H_Rep( cone1 );
    cone2_h := Polymake_H_Rep( cone2 );

    new_ineqs := Concatenation( cone1_h!.inequalities, cone2_h!.inequalities );
    new_equ := Concatenation( cone1_h!.equalities, cone2_h!.equalities );

    return Polymake_ConeFromInequalities( new_ineqs, new_equ );

end );


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

#F qZT( n ) - Z Gate non-terminal
#F Definition: (2^n x 2^n)-matrix  that applies an n-point Z Transform to the target qubits
#F Note:       qZT(1) denotes the matrix[[1, 0], [0, -1]],
#F             qZT(_) is symmetric.
#F
#F qZT( n ) -> an n-qubit Z transform
Class(qZT, TaggedNonTerminal, rec(
  abbrevs   := [ n -> Checked(IsPosInt(n), [n]) ],
  dims      := self >> let(size := 2^self.params[1], [size, size]),
  terminate := self >> Tensor(Replicate(self.params[1], qZ())),
  isReal    := self >> true,
  groups    := self >> [Replicate(self.params[1], 1)],
  rChildren := self >> self.params,
  from_rChildren := (self, rch) >> self.__bases__[1](rch[1]),
  recursive_def := (self, arch) >> self.__bases__[1](self.params[1]),
  SmallRandom := () -> Random([2..5]),
  LargeRandom := () -> Random([6..15]),
  normalizedArithCost := self >> Error("ArithCost not implemented"),
  TType := T_Complex(64)
));

NewRulesFor(qZT, rec(

    # qZT_BinSplit rule
    # qzT_BinSplit qZT_(k) -> (qZT_(k1) tensor I) (I tensor qZT_(k2))
    # k1 + k2 = k
    qZT_BinSplit := rec (
        forTransposition := false,
        minSize          := 2,
        applicable       := (self, nt) >> nt.params[1] > 1,
        children         := nt -> List( [1..nt.params[1] - 1],  i -> [ Tensor(qZT(i), qZT(nt.params[1]-i)).withTags(nt.getTags()) ]  ),
        apply            := (nt, c, cnt) -> c[1],
        switch := true,
    ),

    #F qZT_Base: qZT(1) = qZ() SPL object
    #F Directly represent as an implementable gate
    qZT_Base := rec(
        info             := "qZT_(1) -> qZ()",
        forTransposition := false,
        applicable       := (self, nt) >> nt.params[1]=1,
        apply            := (nt, c, cnt) -> qZ(),
    )

));

ext: gi

max_stars_repo_name: spiral-software/spiral-package-quantum

max_issues_repo_name: spiral-software/spiral-package-quantum

max_forks_repo_name: spiral-software/spiral-package-quantum


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

FFTXGlobals.confGPU := (arg) >> ApplyFunc(arg[1].defaultCUDADeviceConf, Drop(arg,1));
spiral.LocalConfig.fftx := FFTXGlobals;

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:
%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
    /.
        static   tokenKind : number[]=  new Array(128)  ;
         static  __b_init : boolean = %action_type.init_block(%action_type.tokenKind);
        static  init_block(tokenKind : number[]) : boolean
        {
            for (let i = 0; i < tokenKind.length; ++i) {
                tokenKind[i] = 0;
            }
            tokenKind['$'.charCodeAt(0)] = %sym_type.%prefix%DollarSign%suffix%;
            tokenKind['%'.charCodeAt(0)] = %sym_type.%prefix%Percent%suffix%;
            tokenKind['_'.charCodeAt(0)] = %sym_type.%prefix%_%suffix%;

            tokenKind['a'.charCodeAt(0)] = %sym_type.%prefix%a%suffix%;
            tokenKind['b'.charCodeAt(0)] = %sym_type.%prefix%b%suffix%;
            tokenKind['c'.charCodeAt(0)] = %sym_type.%prefix%c%suffix%;
            tokenKind['d'.charCodeAt(0)] = %sym_type.%prefix%d%suffix%;
            tokenKind['e'.charCodeAt(0)] = %sym_type.%prefix%e%suffix%;
            tokenKind['f'.charCodeAt(0)] = %sym_type.%prefix%f%suffix%;
            tokenKind['g'.charCodeAt(0)] = %sym_type.%prefix%g%suffix%;
            tokenKind['h'.charCodeAt(0)] = %sym_type.%prefix%h%suffix%;
            tokenKind['i'.charCodeAt(0)] = %sym_type.%prefix%i%suffix%;
            tokenKind['j'.charCodeAt(0)] = %sym_type.%prefix%j%suffix%;
            tokenKind['k'.charCodeAt(0)] = %sym_type.%prefix%k%suffix%;
            tokenKind['l'.charCodeAt(0)] = %sym_type.%prefix%l%suffix%;
            tokenKind['m'.charCodeAt(0)] = %sym_type.%prefix%m%suffix%;
            tokenKind['n'.charCodeAt(0)] = %sym_type.%prefix%n%suffix%;
            tokenKind['o'.charCodeAt(0)] = %sym_type.%prefix%o%suffix%;
            tokenKind['p'.charCodeAt(0)] = %sym_type.%prefix%p%suffix%;
            tokenKind['q'.charCodeAt(0)] = %sym_type.%prefix%q%suffix%;
            tokenKind['r'.charCodeAt(0)] = %sym_type.%prefix%r%suffix%;
            tokenKind['s'.charCodeAt(0)] = %sym_type.%prefix%s%suffix%;
            tokenKind['t'.charCodeAt(0)] = %sym_type.%prefix%t%suffix%;
            tokenKind['u'.charCodeAt(0)] = %sym_type.%prefix%u%suffix%;
            tokenKind['v'.charCodeAt(0)] = %sym_type.%prefix%v%suffix%;
            tokenKind['w'.charCodeAt(0)] = %sym_type.%prefix%w%suffix%;
            tokenKind['x'.charCodeAt(0)] = %sym_type.%prefix%x%suffix%;
            tokenKind['y'.charCodeAt(0)] = %sym_type.%prefix%y%suffix%;
            tokenKind['z'.charCodeAt(0)] = %sym_type.%prefix%z%suffix%;
            return true;
        }

        public  static    getKind(c :number ):number
        {
            return ((c & 0xFFFFFF80) == 0 /* 0 <= c < 128? */ ? %action_type.tokenKind[c] : 0);
        }
    ./
%End


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

#F --------------------------------------------------------------------------------
RulesFor(DST1, rec(
  #F DST1_Base1: DST1_1 = I_1
  #F
  DST1_Base1 := rec (
      info             := "DST_1 -> I_1",
      forTransposition := false,
      isApplicable     := P -> P[1] = 1,
      rule := (P, C) -> I(1)
  ),

  #F DST1_Base2: DST1_2 = cos(pi/6) * F_2
  #F
  DST1_Base2 := rec (
      info             := "DST_2 -> c * F_2",
      forTransposition := false,
      isApplicable     := P -> P[1] = 2,
      rule := (P, C) -> CosPi(1/6) * F(2)
  ),

  #F DST1_DST3and1: 1984
  #F
  #F   DST1_(n-1) = perm * (DST3_n/2 dirsum DST1_(n/2-1) ^ perm) * (1 tensor DFT_2)^perm
  #F
  #F   Wang:
  #F     Fast Algorithms for the Discrete W Transform and the
  #F     Discrete Fourier Transform.
  #F     IEEE Trans. on ASSP, 1984, pp. 803--814
  #F
  DST1_DST3and1 := rec (
      info          := "DST1_(n-1) --> DST3_n/2, DST1_(n/2-1)",
      isApplicable  := P -> (P[1] + 1) mod 2 = 0 and P[1] >= 2,
      allChildren   := P -> [[ DST3((P[1] + 1)/2), DST1((P[1] - 1)/2) ]],
      rule := (P, C) -> # P is odd
	 DirectSum(C[1], C[2]) ^ OS(P[1], 2) *
	 DirectSum(I(1), Tensor(I((P[1] - 1)/2), F(2))) ^ Z(P[1],P[1]-1) *
	 LIJ(P[1])
         #LIJ(P) *
	 #DirectSum(C[1], C[2] ^ J((P - 1)/2)) *
	 #DirectSum(I(1), Tensor(I((P - 1)/2), F(2))) ^ (Z(P,P-1) * LIJ(P))
  ),

  DST1_toDCT3 := rec( # this rule is probably numerically unstable
      isApplicable  := P -> P[1] >= 3,
      allChildren   := P -> [[ DCT3(P[1]+1) ]],
      rule := (P, C) -> let(n:=P[1], N := n+1,
          RowTensor(n, 1, Mat([[1,-1]])) * C[1] *
          VStack(RowVec(fConst(TReal, n, 0)),
                 Diag(FList(TReal, List([1..n], i -> CosPi(i/(2*N)) / SinPi(i/N))))))
  )
));

#F --------------------------------------------------------------------------------
RulesFor(DST2, rec(
    #F DST2_Base2: DST2_2 = diag(1/Sqrt(2), 1) * F_2
    #F
    DST2_Base2 := rec(
	info := "DST_2 -> F_2",
	isApplicable := P -> P[1] = 2,
	rule := (P, C) -> Diag([Sqrt(1/2), 1]) * F(2)
    ),

    #F DST2_toDCT2: DST2_n =  perm * DCT2_n * diag
    #F
    #F   Wang:
    #F     A Fast Algorithm for the Discrete Sine Transform by the
    #F     Fast Cosine Transform, IEEE Trans. on ASSP, 1982, pp. 814--815
    #F
    DST2_toDCT2 := rec (
	info             := "DST2_n -> DCT2_n",
	isApplicable     := P -> P[1] > 2,
	allChildren := P -> [[ DCT2(P[1]) ]],
	rule := (P, C) ->
	    J(P[1]) * C[1] * Diag(List([0..P[1] - 1], i -> (-1)^i))
    ),

    #F DST2_DST2and4: 1982
    #F
    #F   DST2'_n = perm * (DST4_n/2 dirsum DST2'_n/2^perm) * (1 tensor DFT_2)^perm
    #F
    #F   Wang:
    #F      Fast Algorithms for the Discrete W Transform and the
    #F      Discrete Fourier Transform.
    #F      IEEE Trans. on ASSP, 1984, pp. 803--814
    #F
    DST2_DST2and4 := rec(
	info             := "DST2_n --> DST4_n/2, DST2_n/2",
	isApplicable     := P -> P[1] > 2 and P[1] mod 2 = 0,
	allChildren      := P -> [[ DST4(P[1]/2), DST2(P[1]/2) ]],
	rule := (P, C) ->
	    DirectSum(C[1], C[2]) ^ L(P[1],2) *
	    Tensor(I(Int((P[1] / 2))), F(2)) *
 	    LIJ(P[1])
    )
));

#F --------------------------------------------------------------------------------
RulesFor(DST4, rec(
    #F DST4_Base2: DST4_2 = (1,2) * R_15/8
    #F
    DST4_Base := rec (
	info             := "DST_4 -> R_15/8",
	forTransposition := false,
	isApplicable     := P -> P[1] = 2,
	rule := (P, C) -> J(2) * Rot(15/8)
    ),

    #F DST4_1: 1984
    #F
    #F   DST4_n = perm * DCT4_n * diag
    #F
    #F   Wang:
    #F     Fast Algorithms for the Discrete W Transform and the
    #F     Discrete Fourier Transform.
    #F     IEEE Trans. on ASSP, 1984, pp. 803--814
    #F   Chan:
    #F     Direct Methods for computing discrete sinusoidal transforms
    #F     IEE Proceedings, Vol. 137, 1990, pp. 433--442
    #F
    DST4_toDCT4 := rec (
	info             := "DST4_n --> DCT4_n",
	isApplicable     := P -> P[1] > 2,
	allChildren := P -> [[ DCT4(P[1]) ]],
	rule := (P, C) ->
  	    J(P[1]) * C[1] * Diag(List([0..P[1]-1], i -> (-1)^i))
    )
));


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

################################################################################
##
#F compute_simples(<group>,<int>,<function>) . . . . . compute simples in D^w(G)
##
InstallGlobalFunction( projective_reps, function(G,alpha,n)
  local chi, chi_i, rho_i, Mi, g, x, rho, component_bases, Degrees, algebra_basis_vectors, algebra_basis, twisted_group_alg, dd, lst, i, j, k, h, GSet, T, F;

  F:=CF(n);;

  # cosntruct the twisted group algebra according to the multiplication rule:
  #  g*h=alpha(g,h)gh
  T:=EmptySCTable(Order(G),0);;
  GSet:=AsSet(G);;
  for i in [1..Order(G)] do
    g:=GSet[i];;
    for j in [1..Order(G)] do
      h:=GSet[j];;
      k:=Position(GSet,g*h);;
      # we need to convert alpha(g,h) from an abstract cyclic group element to a
      # complex, so we assign the generator of the cyclic group Z_s to E(s), a
      # primitive s-th root of unity
      #lst:=[E(s)^Position(OrderedA,alpha(g,h)),k];;
      lst:=[alpha(g,h),k];;
      SetEntrySCTable(T,i,j,lst);;
    od;;
  od;;
  twisted_group_alg:=AlgebraByStructureConstants(F,T);;

  # verify that the algebra is associative and then compute its decomposition
  if not(IsAssociative(twisted_group_alg)) then
    Error("algebra is not associative. There is likely an error or your 3-cocycle\n");
  fi;
  dd:=DirectSumDecomposition(twisted_group_alg);

  # compute a basis for the algebra and a basis for each of the decomposition
  # components
  algebra_basis:=Basis(twisted_group_alg);
  algebra_basis_vectors:=BasisVectors(algebra_basis);
  component_bases:=List(dd,Basis);
  Degrees:=List(List(component_bases, Size), Sqrt);# check whether the degrees are squares otherwise extends the field.

  # compute the represenations and characters
  rho:=[];
  chi:=[];
  for Mi in component_bases do
    # compute the representation corresponding to component Mi
    #
    # for each g (basis vector for the twisted group algebra) compute the matrix
    # rho_i(g)_{a,b} according to rho_i(g)(v_a) = Sum_{b}rho_i(g)_{a,b}v_b where
    # {v_{a}} is the basis Mi for the ith component.
    rho_i:=List(algebra_basis_vectors,g->List(Mi,x->Coefficients(Mi,g*x)));
    Append(rho,[rho_i]);
    # for each g, compute chi_i(g)=Tr(rho_i(g))/sqrt(dim of the ideal)
    chi_i:=List(rho_i,x->Trace(x))/Sqrt(Size(Mi));
    Append(chi,[chi_i]);
  od;
  return chi;
end);

#E Projective_Representations.gi . . . . . . . . . . . . . . . . . . . ends here

max_stars_repo_name: pnnl/GroupTheoretical

max_issues_repo_name: pnnl/GroupTheoretical

max_forks_repo_name: pnnl/GroupTheoretical


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

Class(SMPGlobals, rec(
    pthreads := self >> "pthreads",
    OSThreads := self >> Cond(LocalConfig.osinfo.isWindows(), "winthreads", "LinuxThreads"),
    threads := self >> "threads",
    OpenMP := self >> "OpenMP",
    maxThreads := self >> LocalConfig.cpuinfo.cores,
    getOpts := meth(arg)
        local self, opts, optrec, tid;

        self := arg[1];
        opts := CopyFields(SpiralDefaults);
        optrec := rec(api := self.OpenMP(), numproc := self.maxThreads(), parOdd := false);
        if Length(arg) >= 2 then optrec := CopyFields(optrec, arg[2]); fi;

        opts.breakdownRules.GT := Concat([ GT_Base, GT_NthLoop, CopyFields(GT_Par, rec(parEntireLoop := false, splitLoop := true))],
            When(optrec.parOdd, [GT_Par_odd ], []));
        opts.breakdownRules.TTensorI := [ CopyFields(TTensorI_toGT, rec(applicable := (self, t) >> t.hasTags() and ObjId(t.getTags()[1])=AParSMP ))];
        opts.breakdownRules.TTensorInd := [ dsA_base_smp, dsA_smp, L_dsA_L_base_smp, L_dsA_L_smp ];

        tid := When(optrec.api = "OpenMP", threadId(), CopyFields(var("tid", TInt), rec(isParallelLoopIndex := true)));
        opts.tags := [ AParSMP(optrec.numproc, tid) ];

#        if optrec.api = "threads" then opts.subParams := [var("num_threads", TInt), var("tid", TInt)]; fi;
        opts.smp := optrec;

        return opts;
    end
));


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

####################################################
##
## NOTE: Chained unparser!
##

_bitsInt := (bits) -> Sum(List([1..Length(bits)], i -> bits[i]*2^(i-1)));

_sklr_unaligned_store := (bits) -> meth(self, o, i, is )
    local loc0, loc1, offs, p, loffs, mask, l;

    if  IsUnalignedPtrT(o.loc.t) then
        Error("cannot do unaligned store when base pointer is not aligned");
    fi;

    loc0 := vtref(o.exp.t, o.loc, idiv(o.offs, bits));
    loc1 := vtref(o.exp.t, o.loc, idiv(o.offs, bits)+1);
    offs := _unwrap(o.offs) mod bits;
    p    := _unwrap(o.p);
    loffs:= p+offs-bits;
    mask := (n) -> o.exp.t.value(Replicate(n, 1) :: Replicate(bits-n, 0));
    l    := Cond(
                (IsSymbolic(offs) and (p>1)) or (not IsSymbolic(offs) and (offs + p > bits)),
                          # a0 = a0 xor (( a0 xor (val << offs)) & (mask << offs));
                          # a1 = (offs + p <= bits) ? a1 :
                          #          ((a1 >> loffs) << loffs) | ((val & mask) >> bits-offs)
                        [ assign(loc0, bin_xor(loc0, bin_and(bin_shl(mask(p), offs), bin_xor(loc0, bin_shl(o.exp, offs))))),
                          assign(loc1, cond(leq(offs + p, bits), loc1, bin_or(bin_shl(bin_shr(loc1, loffs), loffs), bin_shr(bin_and(o.exp, mask(p)), bits-offs))))],
                p = bits and offs = 0,
                        [ assign(loc0, o.exp) ],
                #else
                        [ assign(loc0, bin_xor(loc0, bin_and(bin_shl(mask(p), offs), bin_xor(loc0, bin_shl(o.exp, offs))))) ]);
    return Print(Blanks(i), self.infix(ESReduce(l, self.opts), Blanks(i)));
end;

Class(SKLRBUnparser, SSEUnparser, rec(
   # none exp must be gone before it reaches unparser
   noneExp := (self, o, i, is) >> self(o.args[1].zero(), i, is),

   utype := (self, t) >> (ObjId(t) = BitVector or (ObjId(t) = TVect and t.t in  [T_Int(1), T_UInt(1), TInt, TUInt]))
                         and t.size in [16, 32, 64],

   BitVector := (self, t, vars, i, is) >> self.T_UInt(T_UInt(t.size), vars, i, is),

   Value := (self, o, i, is) >> When( self.utype(o.t),
        #self(T_UInt(o.t.size).value(_bitsInt(o.v)), i, is),
        Print("0x", HexStringInt(_bitsInt(o.v)), When(o.t.size=64, "LL", "")),
        Inherited(o, i, is)),

   mul := (self, o, i, is) >> Cond( self.utype(o.t), Print("(", self.infix(o.args, ")&("), ")"),
                                    # NOTE: HACK HACK 16x8i mul -> binary and
                                    o.t = TVect(T_Int(8), 16), self( ApplyFunc(bin_and, o.args), i, is ),
                                    Inherited(o, i, is)),
   add := (self, o, i, is) >> When( self.utype(o.t), Print("(", self.infix(o.args, ")^("), ")"), Inherited(o, i, is)),
   sub := ~.add,

   rCyclicShift := (self, o, i, is) >> When( not self.utype(o.t), Inherited(o, i, is),
                      Cond( o.t.size = 16, self.printf("__rotl16($1, $2)", [o.args[1], o.args[2]]),
                            o.t.size = 32, self.printf("__rotl32($1, $2)", [o.args[1], o.args[2]]),
                            o.t.size = 64, self.printf("__rotl64($1, $2)", [o.args[1], o.args[2]]),
                            Error("unsupported size")
                      )),

   sklr_bcast_32x1i := (self, o, i, is) >> self.printf("$1", [tcast(o.t, arith_shr(tcast(T_Int(32), bin_shl(o.args[1], 31-o.args[2])), 31))]),
   sklr_bcast_64x1i := (self, o, i, is) >> self.printf("$1", [tcast(o.t, arith_shr(tcast(T_Int(64), bin_shl(o.args[1], 63-o.args[2])), 63))]),

   sklr_storeu_32x1i := _sklr_unaligned_store(32),
   sklr_storeu_64x1i := _sklr_unaligned_store(64),
));

SKLRBUnparser.bin_shl := (self, o, i, is) >> When( self.utype(o.t) or IsOrdT(o.t),
         self.printf("(($1)<<($2))", [o.args[1], o.args[2]]),
     # else
         Inherited(o, i, is));
SKLRBUnparser.bin_shr := (self, o, i, is) >> Cond(
       self.utype(o.t) or IsUIntT(o.t), self.printf("(($1)>>($2))", [o.args[1], o.args[2]]),
       self.utype(o.t) or IsIntT(o.t),
           Print("((", self.declare(o.t, [], 0, 0), ")(((", self.declare(When(ObjId(o.t)=T_Int, T_UInt(o.t.params[1]), TUInt), [], 0, 0), ")(",
               self(o.args[1], i, is), "))>>(", self(o.args[2], i, is), ")))"),
       Inherited(o, i, is));
SKLRBUnparser.arith_shr := (self, o, i, is) >> Cond(
       IsIntT(o.t), self.printf("(($1)>>($2))", [o.args[1], o.args[2]]),
       self.utype(o.t) or ObjId(o.t)=T_UInt,
           Print("((", self.declare(o.t, [], 0, 0), ")(((", self.declare(T_Int(When(ObjId(o.t)=T_Int,o.t.params[1], o.t.size)), [], 0, 0), ")(",
               self(o.args[1], i, is), "))>>(", self(o.args[2], i, is), ")))"),
       Inherited(o, i, is));
SKLRBUnparser.bin_and := (self, o, i, is) >> When( self.utype(o.t) or o.t in [T_Int(32), T_UInt(32), T_Int(64), T_UInt(64)],
       Print("((", self.infix(o.args, ")&("), "))"),
       Inherited(o, i, is));
SKLRBUnparser.bin_xor := (self, o, i, is) >> When( self.utype(o.t) or o.t in [T_Int(32), T_UInt(32), T_Int(64), T_UInt(64)],
       Print("(", self.infix(o.args, ")^("), ")"),
       Inherited(o, i, is));
SKLRBUnparser.bin_or := (self, o, i, is) >> When( self.utype(o.t) or o.t in [T_Int(32), T_UInt(32), T_Int(64), T_UInt(64)],
       Print("(", self.infix(o.args, ")|("), ")"),
       Inherited(o, i, is));
SKLRBUnparser.tcast := (self, o, i, is) >> When( not self.utype(o.t), Inherited(o, i, is),
       Print("((", self.declare(o.args[1], [], i, is), ")(", self(o.args[2], i, is), "))"));


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

Class(InterpolateSegmentDFT, TaggedNonTerminal, rec(
    abbrevs := [ (numseg, outsize, insize, overlap, dsfunc,usfunc) -> [numseg, outsize, insize, overlap, dsfunc, usfunc] ],

    dims := self >> [self.params[2], self.params[3]],

    terminate := self >> let(
        numseg  := self.params[1],
        rowlength := self.params[2],
        insize  := self.params[3],
        j       := Ind(numseg),
        overlap := self.params[4],
        downsample := self.params[5].at(j),
        upsample   := self.params[6],
        n       := upsample.domain(),
        N       := downsample.range(),
        substrec := i -> rec((j.id):=V(i)),

        kernel := Mat(MatAMat(RowDirectSum(
            overlap,
            List([1..numseg], i ->
                Downsample(SubstVars(Copy(downsample), substrec(i-1))).terminate() *
                DFT(N, 1).terminate() *
                Upsample(upsample).terminate() *
                Scale(1/n, DFT(n, -1).terminate())
            )).toAMat() *
            Scat(fAdd(n*numseg-(numseg-1)*overlap, insize, 0)).toAMat()
        )),
        When(Rows(kernel) = rowlength, kernel, VStack(kernel, O(rowlength-Rows(kernel), Cols(kernel))))
    ),
    isReal    := False,
    normalizedArithCost := self >> let(numseg := self.params[1], n:=self.params[6].domain(), N := self.params[6].range(),
        numseg * (n + IntDouble(5 * n * d_log(n) / d_log(2)) + IntDouble(5 * N * d_log(N) / d_log(2)))),
    TType := T_Complex(T_Real(64)),

    HashId := self >> let(h := [ self.params[1], self.params[2], self.params[3], self.params[4], self.params[5].domain(), self.params[5].range()],
        When(IsBound(self.tags), Concatenation(h, self.tags), h)),
    doNotMeasure := true
));

_sumSegDims := function(numsegs, seglen)
    local l, lfact, lsum, divsum, ssum;

    l := List([0..numsegs-1], i->seglen.at(i).domain());

    # check to pull out the mul
    if ForAll(l, i->ObjId(i) = mul and IsValue(i.args[1]) and i.args[1].v = l[1].args[1].v) then
        lfact := l[1].args[1];
        lsum := List(l, i -> i.args[2]);
        # we may need to create a table, but not sure yet...
        divsum := ApplyFunc(add, lsum);
        ssum := mul(lfact, divsum);
    else
        ssum := ApplyFunc(add, l);
    fi;

    return ssum;
end;

NewRulesFor(InterpolateSegmentDFT, rec(
    InterpolateSegmentDFT_base := rec(
        switch := false,
        applicable := (self, nt) >> not nt.hasTags() and ObjId(nt.params[6]) = fZeroPadMiddle,
        children := nt -> let(numseg := nt.params[1], j := Ind(numseg), downsample := nt.params[5].at(j), upsample := nt.params[6],
                              insize := nt.params[3], overlap := nt.params[4], n:=upsample.domain(), N := downsample.range(),
                              us := N/n, blk := n/2,
                              inp := n + (numseg-1)*(n-overlap), over := inp - insize, l := n - over,
            [[
                Downsample(downsample),
                DFT(N, 1),
                Upsample(upsample),
                DFT(n, -1),
                Upsample(fAdd(n, l, 0)),
                InfoNt(j)
            ]]),
        apply := (nt, C, cnt) -> let(numseg := nt.params[1], _j := cnt[6].params[1], j := Ind(_j.range-1), downsample := nt.params[5].at(j), upsample := nt.params[6],
                              outsize := nt.params[2], insize := nt.params[2], overlap := nt.params[4], n:=upsample.domain(), N := downsample.range(),
                              ids_rows := _sumSegDims(numseg-1, nt.params[5]),
                RowDirectSum(overlap,[
                    IRowDirSum(j, numseg-1, overlap,
                        SubstVars(Copy(C[1]), rec((_j.id) := j)) * C[2] * C[3] * Scale(1/n, C[4])).overrideDims([ids_rows, insize - C[4].dims()[2] + overlap]),
                    SubstVars(Copy(C[1]), rec((_j.id) := V(numseg-1))) * C[2] * C[3] * Scale(1/n, C[4] * C[5])
                ]).overrideDims([nt.params[2], nt.params[3]])
       )
    ),
    InterpolateSegmentDFT_PrunedDFT := rec(
        applicable := (self, nt) >> not nt.hasTags() and ObjId(nt.params[6]) = fZeroPadMiddle,
        children := nt -> let(numseg := nt.params[1], j := Ind(numseg), downsample := nt.params[5].at(j), upsample := nt.params[6],
                              insize := nt.params[3], overlap := nt.params[4], n:=upsample.domain(), N := downsample.range(),
                              us := N/n, blk := n/2,
                              inp := n + (numseg-1)*(n-overlap), over := inp - insize, l := n - over, g := Gcd(n, l),
            [[
                Downsample(downsample),
                PrunedDFT(N, 1, blk, [0, 2*us-1]),
                DFT(n, -1),
                PrunedDFT(n, -1, g, [0..l/g-1]),
                InfoNt(j)
            ]]),
        apply := (nt, C, cnt) -> let(numseg := nt.params[1], _j := cnt[5].params[1], j := Ind(_j.range-1), downsample := nt.params[5].at(j), upsample := nt.params[6],
                              outsize := nt.params[2], insize := nt.params[3], overlap := nt.params[4], n := upsample.domain(), N := downsample.range(),
                              ids_rows := _sumSegDims(numseg-1, nt.params[5]),
                RowDirectSum(overlap,[
                    IRowDirSum(j, numseg-1, overlap,
                        SubstVars(Copy(C[1]), rec((_j.id) := j)) * C[2] * Diag(fConst(n, 1/n)) * C[3]).overrideDims([ids_rows, insize - C[4].dims()[2] + overlap]),
                    SubstVars(Copy(C[1]), rec((_j.id) := V(numseg-1))) * C[2] * Diag(fConst(n, 1/n)) * C[4]
                ]).overrideDims([nt.params[2], nt.params[3]])
       )
    )
));


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

Class(CudaUnparser, CUnparser, rec(

    simt_loop := meth(self, o, i, is)
        local v, lo, hi, l_cond, r_cond, simt_idx, simt_rng, body;

        v := o.var; lo := o.range[1]; hi := Last(o.range);
        simt_idx := o.simt_idx;
        simt_rng := o.simt_idx.get_rng();
        body := Copy(o.cmd);
        l_cond := simt_rng[1]>0; r_cond := simt_rng[2] < simt_idx.size()-1;

        if (hi-lo) < simt_rng[2]-simt_rng[1] then
            body := SubstVars(body, rec((v.id) := simt_idx + lo - simt_rng[1]));
            Print(Blanks(i), "if( ");
            if l_cond then
                Print(self(simt_idx, 0, 0), " >= ", self(simt_rng[1], 0, 0), " && ");
            fi;
            Print(self(simt_idx, 0, 0), " <= ", self(simt_rng[1]+hi-lo, 0, 0), " ) {\n",
                self(body, i+is, is),
                Blanks(i), "}\n");
        elif (hi-lo) = simt_rng[2]-simt_rng[1] then
            body := SubstVars(body, rec((v.id) := simt_idx + lo - simt_rng[1]));
            if l_cond or r_cond then
                Print(Blanks(i), "if( ");
                if l_cond then
                    Print(self(simt_idx, 0, 0), " >= ", self(simt_rng[1], 0, 0));
                fi;
                if l_cond and r_cond then
                    Print(" && ");
                fi;
                if r_cond then
                    Print(self(simt_idx, 0, 0), " <= ", self(simt_rng[2], 0, 0));
                fi;
                Print(" ) {\n");
            fi;
            self(body, When(l_cond or r_cond, i+is, i), is);
            if l_cond or r_cond then
                Print(Blanks(i), "}\n");
            fi;
        else
            body := SubstVars(body, rec((v.id) := simt_idx + v - simt_rng[1]));
            Print(When(IsBound(self.opts.looppragma), self.opts.looppragma(o,i,is)),
                Blanks(i), "for(int ", v, " = ", lo, "; ");
            if l_cond then
                Print(self(simt_idx + v-simt_rng[1], 0, 0), " >= ", lo, " && ");
            fi;
            Print(self(simt_idx + v-simt_rng[1], 0, 0), " <= ", hi, "; ",
                v, "+=", self(simt_rng[2]-simt_rng[1]+1, 0, 0), ") {\n",
                self(body, i+is, is),
                Blanks(i), "}\n");
        fi;

    end,

    Dim3 := meth(self,t,vars,i,is)
        local v, var_list;
        Print("dim3 ");
        var_list := List(vars, v-> v.id::"("::StringJoin(", ", List([v.x.value, v.y.value, v.z.value], v -> When(IsValue(v), v.v, v) ) )::")");
        self.infix(var_list, ", ");
    end,

    TArray := meth(self,t,vars,i,is)
        local dims, elt, v, ptype, vsize, var_specs;
        if Length(vars) > 1 then DoForAll(vars, v->Print(self.TArray(t, [v], i, is), "; "));
        elif Length(vars) = 0 then
            Print(self.declare(t.t, [], i, is), " *");
        else
            # at this point Length(vars)=1
            v := When(IsList(vars), vars[1], vars);
            var_specs := When(IsBound(v.decl_specs), v.decl_specs, []);
            dims := []; elt := t;
            while IsArray(elt) do Add(dims, elt.size); elt := elt.t; od;

            #FIXME/HACK: Ignoring twiddles by looking for "D" in .id
            # Better way: look for func.name="init" in parent/context
            if IsBound(self.opts.useMemoryArena) and self.opts.useMemoryArena and v.id[1] <> 'D' then
                #FIXME: Arena currently doesn't handle multiple dims.
                #FIXME: Slightly ugly hack to get this to be a pointer

                # To handle vectors. Arena is declared only for scalars. For
                # vectors, we must manually scale the allocation by vector length.
                vsize := 1; if ObjId(elt)=TVect then vsize := elt.size; fi;
                ptype := Concatenation(elt.name, "Pointer");
                self.(ptype)(elt, [v], i, is);
                Print(" =  &(ARENA[ (arenalevel-=",self(dims[1]*vsize,i,is),") ])");
            else
                if var_specs <> [] then
                    Print(self.infix(var_specs, " "), " ");
                fi;
                if IsBound(v.decl_as_ptr) and v.decl_as_ptr then
                    Print(self.declare(t.t, [], i, is), "*", v.id);
                else
                    self.(elt.name)(elt, [v], i, is);
                    DoForAll(dims, d->Print("[",self(d,i,is),"]"));
                fi;
            fi;

        fi;
    end,

    T_CudaEvent := (self, t, vars, i, is) >>
        Print("cudaEvent_t ", self.infix(vars, ", ", i + is)),

    cu_event_create := (self, o, i, is) >> Print(Blanks(i), "cudaEventCreate(", self(o.args[1], 0, 0),
                            When(Length(o.args)>1 and not o.args[2].v, ")", ");\n")),

    cu_event_record := (self, o, i, is) >> Print(Blanks(i), "cudaEventRecord(", self(o.args[1], 0, 0),
                            When(Length(o.args)>1 and not o.args[2].v, ")", ");\n")),

    cu_event_destroy := (self, o, i, is) >> Print(Blanks(i), "cudaEventDestroy(", self(o.args[1], 0, 0),
                            When(Length(o.args)>1 and not o.args[2].v, ")", ");\n")),

    cu_event_elapsed_time := (self, o, i, is) >> Print(Blanks(i), "cudaEventElapsedTime",
                                self.pinfix(o.args{[1..3]}, ", "),
                                When(Length(o.args)>3 and not o.args[4].v, "", ";\n")),

    cu_event_synchronize := (self, o, i, is) >> Print(Blanks(i), "cudaEventSynchronize(", self(o.args[1], 0, 0),
                            When(Length(o.args)>1 and not o.args[2].v, ")", ");\n")),

    cu_device_synchronize := (self, o, i, is) >> Print(Blanks(i), "cudaDeviceSynchronize(",
                            When(Length(o.args)>0 and not o.args[1].v, ")", ");\n")),

    cprintf := (self, o, i, is) >> Print(Blanks(i), "printf(\"", o.args[1].v, "\", ",
                                            self.infix(o.args{[2..Length(o.args)]}, ", "), ");\n"),

    func := (self, o, i, is) >> let(
        parameters:=Flat(o.params),
        id := Cond(o.id="transform" and IsBound(self.opts.subName),
                     self.opts.subName,
                   o.id="init"      and IsBound(self.opts.subName),
                     Concat("init_",self.opts.subName),
                   o.id="destroy"      and IsBound(self.opts.subName),
                     Concat("destroy_",self.opts.subName),
                   o.id),
		When ((IsBound(self.opts.wrapCFuncs) and self.opts.wrapCFuncs), Print("\nextern \"C\" {")),
        Print("\n", Blanks(i),
            self.opts.funcModifier, self.declare(o.ret, var(id, o.ret), i, is), "(",
            DoForAllButLast(parameters, p->Print(self.declare(p.t, p,i,is), ", ")),
            When(Length(parameters)>0, self.declare(Last(parameters).t, Last(parameters),i,is), ""), ") ",
            "{\n",
            When(IsBound(self.opts.postalign), DoForAll(parameters, p->self.opts.postalign(p,i+is,is))),
            self(o.cmd, i+is, is),
            Blanks(i),
            "}\n"),
		When ((IsBound(self.opts.wrapCFuncs) and self.opts.wrapCFuncs), Print("}\n"))),

    specifiers_func := (self, o, i, is) >> let(
        parameters:=Flat(o.params),
        id := o.id,
        Print("\n", Blanks(i),
            self.opts.funcModifier, self.infix(o.decl_specs, " "), " ", self.declare(o.ret, var(id, o.ret), i, is), "(",
            DoForAllButLast(parameters, p->Print(self.declare(p.t, p,i,is), ", ")),
            When(Length(parameters)>0, self.declare(Last(parameters).t, Last(parameters),i,is), ""), ") ",
            "{\n",
            When(IsBound(self.opts.postalign), DoForAll(parameters, p->self.opts.postalign(p,i+is,is))),
            self(o.cmd, i+is, is),
            Blanks(i),
            "}\n")),

    cu_call := (self, o, i, is) >>
        Print(Blanks(i), o.func,
                "<<<",
                self.infix([o.dim_grid, o.dim_block], ", "),
                ">>>",
                self.pinfix(o.args, ", "), ";\n"),

    cu_allocate := (self, o, i, is) >> Print(Blanks(i),
        "cudaMalloc(&",
        self(o.loc, 0, 0), ", ", self(o.size, 0, 0), "*sizeof(", self.declare(o.type,[], 0, 0), "));\n"),

    cu_allocate_managed := (self, o, i, is) >> Print(Blanks(i),
        "cudaMallocManaged(&",
        self(o.loc, 0, 0), ", ", self(o.exp.size, 0, 0), "*sizeof(", self.declare(o.exp.t,[], 0, 0), "));\n"),

    cu_memcpy := (self, o, i, is) >> Print(Blanks(i),
        "cudaMemcpy(", self(o.loc, 0, 0), ", ", self(o.exp, 0, 0), ", ", self(o.size, 0, 0), "*sizeof(", self.declare(o.loc.t.t,[], 0, 0),
            "), ", o.kind, ");\n"),

    cu_memcpy_to_sym := (self, o, i, is) >> Print(Blanks(i),
        "cudaMemcpyToSymbol(", self(o.loc, 0, 0), ", ", self(o.exp, 0, 0), ", ", self(o.size, 0, 0), "*sizeof(", self.declare(o.loc.t.t,[],0,0), ") );\n"),

    cu_check_errors := (self, o, i, is) >> Print(Blanks(i), "checkCudaErrors( ", self(o.args[1], 0, 0), " );\n"),

    cu_free := (self, o, i, is) >> Print(Blanks(i), "cudaFree(", self(o.args[1], 0, 0), ");\n"),

#    cospi := (self,o,i,is) >> Print("_"::o.name, self.pinfix(o.args, ", ")),
#    sinpi := (self,o,i,is) >> Print("_"::o.name, self.pinfix(o.args, ", ")),

    simtThreadIdxX := (self, o, i, is) >> Print("threadIdx.x"),
    simtThreadIdxY := (self, o, i, is) >> Print("threadIdx.y"),
    simtThreadIdxZ := (self, o, i, is) >> Print("threadIdx.z"),

    simtBlockIdxX := (self, o, i, is) >> Print("blockIdx.x"),
    simtBlockIdxY := (self, o, i, is) >> Print("blockIdx.y"),
    simtBlockIdxZ := (self, o, i, is) >> Print("blockIdx.z"),

    simt_syncgrid := (self, o, i, is) >> Error("Grid-level sync not supported."),
    simt_syncblock := (self, o, i, is) >> Print(Blanks(i), "__syncthreads();\n"),
    simt_synccluster := (self, o, i, is) >> Print(Blanks(i), "__syncwarp();\n")

    ));

CudaDefaults := CopyFields(SpiralDefaults,
                                rec(
                                    globalUnrolling := 10,
                                    useDeref := false,
                                    generateInitFunc := false,
                                    includes := [],
                                    arrayBufModifier := "",
                                    arrayDataModifier := "",

                                    mempool := true,
                                    use_shmem := true,
                                    gpu_timing := true,
                                    sumsgen := SIMTSumsGen,
                                    codegen :=  CudaCodegen,
                                    unparser := CudaUnparser,
                                    cuda_version := "10.1",
                                )
                    );

TitanVDefaults := CopyFields(CudaDefaults,
                                rec(
                                    codename := "GV100-400-A1",
                                    max_l1_size := 128*1024,
                                    max_shmem_size := 96*1024,
                                    l2_size := 4.5*2^20,
                                    devmem_size := 12*2^30
                                )
                    );

max_stars_repo_name: mikefranusich/spiral-package-simt

max_issues_repo_name: mikefranusich/spiral-package-simt

max_forks_repo_name: mikefranusich/spiral-package-simt


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

# PR(<dom>, <m>, <n>, <j>) - scatter function for our older PRDFT agorithm,
#                            still used in PRDFT rule for Circulant
Class(PR, FuncClass, rec(
    def := (dom,m,n,j) -> rec(N := approx.FloorRat(m*n/2)+1, n := dom),
    lambda := self >> let(
        m := self.params[2], n := self.params[3], j := self.params[4],
        i := Ind(self.params[1]),
        Lambda(i, cond(leq(i, approx.CeilingRat(m/2)-1), i*n + j, m*n - i*n - j)))
));

#####################
# Matrices
#####################

Class(BotHalf, Sym, rec(
    def := (m, conj) -> let(mc := idiv(m+1, 2), mf := idiv(m, 2),
    DirectSum(I(2*mc), Tensor(I(mf), conj)))
));
Class(TopHalf, Sym, rec(
    def := (m, conj) -> let(mc := idiv(m+1, 2), mf := idiv(m, 2),
    DirectSum(Tensor(I(mc), conj), I(2*mf)))
));

Class(BHD, DiagFunc, rec(
    def := (m, e0, e1) -> Checked(IsPosIntSym(m), rec(size := 2*m)),
    domain := self >> 2*self.params[1],
    range := self >> TReal,
    lambda := self >> let(
	m := self.params[1], mc := idiv(m+1, 2), mf := idiv(m, 2),
	e0:=self.params[2], e1:=self.params[3], i := Ind(2*m),
	Lambda(i, cond(leq(i, 2*mc-1), 1.0, cond(neq(imod(i,2),0), e1, e0))))
));

Class(BHN, DiagFunc, rec(
    def := (m) -> Checked(IsPosIntSym(m), rec(size := m)),
    domain := self >> self.params[1],
    range := self >> TReal,
    lambda := self >> let(
    m := self.params[1], mc := idiv(m+1, 2), i := Ind(m),
    Lambda(i, cond(leq(i, mc-1), 1, -1)))
));

Class(HTwid, Sym, rec(
    def := (N, n, k, a, b, i, hconj_mat) -> let(
    na := Numerator(a),   nb := Numerator(b),
    da := Denominator(a), db := Denominator(b),
    w := Global.Conjugate(E(N*da*db)^k),
    DirectSum(List([0..n-1],
        j -> hconj_mat * Mat(RealMatComplexMat([[w ^ ((i*da+na)*(j*db+nb))]])) )))
));

X13 := (n,blk) -> DirectSum(I(1), Tensor(I(n-1), blk)) ^ (Z(2*n-1,-1)*LIJ(2*n-1));
SymSplit1 := n -> DirectSum(I(n+1), J(n-1)) * DirectSum(I(1), X13(n,F(2)));
SymSplit3 := n -> IJ(2*n,n)   * DirectSum(I(1), X13(n,J(2)*F(2)));

X24 := (n,blk) -> Tensor(blk, I(n)) * IJ(2*n,n);
SymSplit2 := n -> X24(n,F(2));
SymSplit4 := n -> X24(n,J(2)*F(2));

rpdiag := p -> Tensor(pdiag(Int(p/2)+1), I(2));
rperm_ev := p -> Tensor(J((p+1)/2), Diag(1,-1));
rperm_od := p -> Tensor(J((p+1)/2), Diag(-1,1));

# YSV: Many rules below turned out to be invalid for odd stride in H, so I added guards
#
RewriteRules(RulesFuncSimp, rec(
 BHD_Stride_fId := ARule(fCompose,
     [ [BHD, @(1),@(2),@(3)], [fTensor, [H,@,@(4),@,@.cond(e->IsEvenInt(e) and e>1)], [fId, 2]] ],
     e -> [BHD(@(4).val, @(2).val, @(3).val)] ),

 BHD_HH_fId := ARule(fCompose,
     [ [BHD, @(1),@(2),@(3)], [fTensor, [HH,@,@(4),_0,[ListClass, @.cond(e->IsEvenInt(e) and e>1), _1]], [fId, 2]] ],
     e -> [BHD(@(4).val, @(2).val, @(3).val)] ),

 BHD_Stride_H := ARule(fCompose,
     [ [BHD, @(1),@(2),@(3)], [fTensor, [H,@,@(4),@,@.cond(e->IsEvenInt(e) and e>1)], [H,2,2,0,1]] ],
     e -> [BHD(@(4).val, @(2).val, @(3).val)] ),

 BHN_H := ARule(fCompose,
     [ BHN, [@(1,H),@,@,@, @.cond(e->IsEvenInt(e) and e>1)] ],
     e -> [BHN(@(1).val.domain())] ),

# ******* TOGGLE THIS RULE
 # NOTE: IsSymbolic in the condition is dangerous! this rule is not valid if stride is odd.
 #        which is precisely what happens with ASPF_CT1Odd_RFT rule, but this rule is still needed for DCTs
# BHN_HH := ARule(fCompose,
#     [ BHN, [@(1,HH),@,@,_0,[ListClass, @.cond(e->IsSymbolic(e) or (IsEvenInt(e) and e>1)), _1]] ], # vector stride = 1 means stride <> 1
#     e -> [BHN(@(1).val.domain())] ),
##############

# NOTE: breaks for odd sizes -> guard
 BHN_fTensor := ARule(fCompose,
     [ BHN, [@(1, fTensor), @(2,fId,e->IsEvenInt(e.range())), @(3,fBase,e->IsEvenInt(e.range()))] ],
     e -> [BHN(@(1).val.domain())] ),

 BHD_fTensor := ARule(fCompose,
     [ @(0,BHD), [@(1, fTensor), @(2,fId,e->IsEvenInt(e.range())), @(3,fBase,e->IsEvenInt(e.range())), [fId, 2]] ],
     e -> [BHD(@(2).val.params[1], @(0).val.params[2], @(0).val.params[3])] ),

 ## BH o H -> BH
 BH_H := ARule(fCompose, [ [@(1,BH), @N, @r, @n, @b, @s], [H, @n, @m, @bb, @(2).cond(e->e>1 and IsEvenInt(e))] ],
     (e, cx) -> [ BH(@N.val, @r.val, @m.val,
                 memo(cx, "g", @b.val + @s.val * @bb.val), # new base
             @s.val * @(2).val) ]),                     # new stride

 ## H o BH -> BH
 H_BH := ARule(fCompose, [ [H, @N, @n, @b, @(2).cond(e->e>1)], [@(1,BH), @n, @r, @m, @bb, @ss] ],
     (e, cx) -> [ When(@b.val<>0, H(@N.val,@N.val,@b.val, 1), fId(@N.val)),
              BH(@N.val, @r.val * @(2).val, @m.val,
                 memo(cx, "g", @(2).val * @bb.val), # new base
             @(2).val * @ss.val) ]),            # new stride

));


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

update_type := function(vars, newType)
  vars.t := newType;
end;

Unique := (lst) ->
       FoldL(lst, (b, a)->Concat( b, When(not a in b, [a], [])), []);

SubstList:= function(expr, list_match, list_replace)
   local new_expr, index;
   new_expr := expr;
   for index in [1..Length(list_match)] do
     new_expr := SubstBottomUp(new_expr, @(2).cond(dc->dc=list_match[index]), x->list_replace[index]);
   od;
   return new_expr;
end;

Class(RulesStateHCOL, RuleSet, rec(inType := "SigmaSPL", outType := "SigmaSPL"));
RewriteRules(RulesStateHCOL, rec(
   Constant_PointWise := ARule(OLCompose, [@(1, PointWise, e->Length(Collect(e.op.expr, e.op.vars[1]))=0), @(2)], e->[@(1).val])
));

Class(RulesTerminateReductionHCOL, RuleSet, rec(inType := "SigmaSPL", outType := "SigmaSPL"));
RewriteRules(RulesTerminateReductionHCOL, rec(
     Reduction_GathH := ARule(OLCompose, [@(1, Reduction), @(2, GathH)],
       e->[ let(o:= @(1).val, i:= Ind(o.N),
        ISumReduction(i, o.N, o.op, o.idval, o.isSaturated, eT(@(2).val.N, add(i*@(2).val.stride, @(2).val.base)))) ]),

    Reduction_terminate := Rule(@(1, Reduction), e->let(o:= @(1).val, i := Ind(o.N),
        ISumReduction(i, o.N, o.op, o.idval, o.isSaturated, eT(o.N, i)))),

	Scat1Union_Gath1 := ARule(OLCompose, [@(1, ScatHUnion, e1->e1.n=1), @(2, GathH, e2->e2.n=1)],
		e->[ let(
		Print("Scat1Union_Gath1: ", @(1).val.N, " ", @(1).val.base, " ", @(2).val.N, " ", @(2).val.base,"\n"),
		OLCompose(
		  eUnion(@(1).val.N, @(1).val.base),
	      eT(@(2).val.N, @(2).val.base)
		)
	)]),

    ScatHUnion_GathH := ARule(OLCompose, [@(1, ScatHUnion), @(2, GathH)],
        e->[let( i:= Ind(@(1).val.n),
	  ISumUnion(i, @(1).val.n,
	    OLCompose(
	      eUnion(@(1).val.n, @(1).val.base+@(1).val.stride*i),
	      eT(@(2).val.n, @(2).val.base+@(2).val.stride*i)
	    )
	  )
	)]),

    GathH0_terminate := Rule(@(1, GathH, e->e.n=0), e->eT(0, -1)),
    GathH1_terminate := Rule(@(1, GathH, e->e.n=1), e->eT(@(1).val.N, @(1).val.base)),
    GathHN_terminate := Rule(@(1, GathH, e->e.n>1),
        e-> let(i := Ind(@(1).val.n),
	    	  ISumUnion(i, @(1).val.n,
		    OLCompose(
			eUnion(@(1).val.n, i),
			eT(@(1).val.N, @(1).val.base+@(1).val.stride*i))))),

    ScatHUnion0_terminate := Rule(@(1, ScatHUnion, e->e.n=0), e->eUnion(0,-1)),   #eUnion(@(1).val.N, @(1).val.base)),
    ScatHUnion1_terminate := Rule(@(1, ScatHUnion, e->e.n<=1), e->eUnion(@(1).val.N, @(1).val.base)),
    ScatHUnionN_terminate := Rule(@(1, ScatHUnion, e->e.n>1),
        e-> let(i := Ind(@(1).val.n),
	       ISumUnion(i, @(1).val.n,
 	         OLCompose(
			eUnion(@(1).val.N, @(1).val.base+@(1).val.stride*i),
			eT(@(1).val.n, i)
	         )
               ))
    ) ));

Class(RulesSumsHCOLv2a, RuleSet, rec(inType := "SigmaSPL", outType := "SigmaSPL"));
RewriteRules(RulesSumsHCOLv2a, rec(
    OLCompose_Assoc := ARule(OLCompose, [ @(1,OLCompose) ],  e -> @(1).val.children() ),
    OLCompose_PointWise_PointWise := ARule(OLCompose, [ @(1, PointWise), @(2, PointWise) ],
        e -> [ PointWise(@(1).val.N, Lambda(@(2).val.op.vars, SubstVars(Copy(@(1).val.op.expr),
         rec((@(1).val.op.vars[2].id) := @(2).val.op.vars[2], (@(1).val.op.vars[1].id) := @(2).val.op.expr)))) ]),
    GathH_GathH := ARule(OLCompose, [@(1, GathH), @(2, GathH)],
        e -> [GathH(@(2).val.N, @(1).val.n, @(1).val.base+@(2).val.base, @(1).val.stride*@(2).val.stride)]),
    Reduction_ScatHUnion := ARule(OLCompose, [ @(1, Reduction), @(2, ScatHUnion, e->e.n=1) ], e->[]),
    PointWise_BinOp := ARule(OLCompose, [@(1, PointWise, e->e.N=1), @(2, BinOp, e->e.N=1)],
        e -> [BinOp(1,Lambda(@(2).val.op.vars, @(1).val.op.at(@(2).val.op.expr, V(0))))]),
    PointWise_ISumUnion :=  ARule(OLCompose, [ @(1, PointWise), @(2, ISumUnion) ],
        e -> [ CopyFields(@(2).val, rec(
            _children :=  List(@(2).val._children, c -> OLCompose(@(1).val, c)),
            dimensions := [Rows(@(1).val), Cols(@(2).val)] )) ]),
    PointWise_ScatHUnion := ARule(OLCompose, [@(1,PointWise), @(2, ScatHUnion)],
        e -> let(i := Ind(@(2).val.dims()[2]), [@(2).val, PointWise(@(2).val.dims()[2],
            Lambda([@(1).val.op.vars[1], i], SubstVars(@(1).val.op.expr, rec((@(1).val.op.vars[2].id) := i*@(2).val.stride+@(2).val.base))))])),
    ISumXXX_YYY := ARule(OLCompose, [ @(1, [ISumUnion, ISumReduction]), @(2, [GathH, PointWise, Induction, eT]) ],
        e -> [ CopyFields(@(1).val, rec(
            _children :=  List(@(1).val._children, c -> OLCompose(c, @(2).val)),
            dimensions := [Rows(@(1).val), Cols(@(2).val)] )) ]),
    SUMUnion_GathH := ARule(OLCompose, [@(1, SUMUnion), @(2, GathH)],
       e ->[SUMUnion(List(@(1).val._children, c->OLCompose(c, @(2).val)))] ),
    ScatHUnion_SUMUnion := ARule(OLCompose, [ @(1, ScatHUnion), @(2, SUMUnion) ],
        e -> [SUMUnion(List(@(2).val._children, c->OLCompose(@(1).val,c)))] ),
    ScatHUnion_ISumUnion := ARule(OLCompose, [@(1, ScatHUnion), @(2, ISumUnion)],
        e -> [ISumUnion(@(2).val.var, @(2).val.domain, OLCompose(@(1).val, @(2).val._children[1]))]),
    ScatHUnion_ScatHUnion := ARule(OLCompose, [@(1, ScatHUnion), @(2, ScatHUnion)],
        e -> [ScatHUnion(@1.val.N, @(2).val.n, @(1).val.base+@(2).val.base, @(2).val.stride)] ),
    SUMUnion_Assoc := ARule(SUMUnion, [@(1, SUMUnion) ], e->@(1).val.children())
));

Class(RulesSumsHCOLv2b, RuleSet, rec(inType := "SigmaSPL", outType := "SigmaSPL"));
RewriteRules(RulesSumsHCOLv2b, rec(
    OLCompose_Assoc := ARule(OLCompose, [ @(1,OLCompose) ],  e -> @(1).val.children() ),
    OLCompose_PointWise_PointWise := ARule(OLCompose, [ @(1, PointWise), @(2, PointWise) ],
        e -> [ PointWise(@(1).val.N, Lambda(@(2).val.op.vars, SubstVars(Copy(@(1).val.op.expr),
         rec((@(1).val.op.vars[2].id) := @(2).val.op.vars[2], (@(1).val.op.vars[1].id) := @(2).val.op.expr)))) ]),
    GathH_GathH := ARule(OLCompose, [@(1, GathH), @(2, GathH)],
        e -> [GathH(@(2).val.N, @(1).val.n, @(1).val.base+@(2).val.base, @(1).val.stride*@(2).val.stride)]),
    ISumReduction_PointWise := ARule(OLCompose, [ @(1, ISumReduction), @(2, PointWise) ],
        e -> [ ISumReduction(@(1).val.var, @(1).val.domain, @(1).val.op, @(1).val.idval, @(1).val.isSaturated,
                OLCompose(@(1).val._children[1], @(2).val))]),

    eT_Pointwise := ARule(OLCompose, [@(1, eT, e->IsVar(e.base)), @(2, PointWise)],
         e -> let(i := Ind(1),
	     [PointWise(1, Lambda(
                SubstVars(Copy(@(2).val.op.vars), rec((@(2).val.op.vars[2].id):=i)),
                SubstVars(Copy(@(2).val.op.expr), rec((@(2).val.op.vars[2].id):=@(1).val.base)))), @(1).val])),
    ISumXXX_YYY := ARule(OLCompose, [ @(1, [ISumUnion, ISumReduction]), @(2, [GathH, PointWise, Induction, eT]) ],
        e -> [ CopyFields(@(1).val, rec(
            _children :=  List(@(1).val._children, c -> OLCompose(c, @(2).val)),
            dimensions := [Rows(@(1).val), Cols(@(2).val)] )) ]),
    eT_ISumUnion := ARule(OLCompose, [ @(1, eT), @(2, ISumUnion, e->ObjId(e._children[1]._children[1])=eUnion)],
        e -> [ Drop(SubstVars(Copy(@(2).val._children[1]._children), rec((@(2).val.var.id) := (@(1).val.base)) ), 1)]),
    eT_Induction := ARule(OLCompose, [@(1, eT), @(2, Induction)],
        e -> [Inductor(@(2).val.N, @(1).val.base, @(2).val.op, @(2).val.initval)]),
  PointWise_ISumUnion :=  ARule(OLCompose, [ @(1, PointWise), @(2, ISumUnion) ],
        e -> [ CopyFields(@(2).val, rec(
            _children :=  List(@(2).val._children, c -> OLCompose(@(1).val, c)),
            dimensions := [Rows(@(1).val), Cols(@(2).val)] )) ]),

   PointWise := ARule(OLCompose, [@(1, PointWise, e->Length(Collect(e.op.expr, e.op.vars[1]))=0), ...],
       e -> [@(1).val])
));

Class(RulesSumsHCOL, RuleSet, rec(inType := "SigmaSPL", outType := "SigmaSPL"));
RewriteRules(RulesSumsHCOL, rec(
    OLCompose_Assoc := ARule(OLCompose, [ @(1,OLCompose) ],  e -> @(1).val.children() ),
    OLCompose_PointWise_PointWise := ARule(OLCompose, [ @(1, PointWise), @(2, PointWise) ],
        e -> [ PointWise(@(1).val.N, Lambda(@(2).val.op.vars, SubstVars(Copy(@(1).val.op.expr),
         rec((@(1).val.op.vars[2].id) := @(2).val.op.vars[2], (@(1).val.op.vars[1].id) := @(2).val.op.expr)))) ]),
    PointWise_ISumUnion :=  ARule(OLCompose, [ @(1, PointWise), @(2, ISumUnion) ],
        e -> [ CopyFields(@(2).val, rec(
            _children :=  List(@(2).val._children, c -> OLCompose(@(1).val, c)),
            dimensions := [Rows(@(1).val), Cols(@(2).val)] )) ]),
    PointWise_ScatHUnion := ARule(OLCompose, [@(1,PointWise), @(2, ScatHUnion)],
        e -> let(i := Ind(@(2).val.dims()[2]), [@(2).val, PointWise(@(2).val.dims()[2],
            Lambda([@(1).val.op.vars[1], i], SubstVars(@(1).val.op.expr, rec((@(1).val.op.vars[2].id) := i*@(2).val.stride+@(2).val.base))))])),
    Reduction_ISumReduction :=  ARule(OLCompose, [ @(1, Reduction), @(2, ISumUnion) ],
        e -> [ ISumReduction(@(2).val.var, @(2).val.domain, @(1).val.op, @(1).val.idval, @(1).val.isSaturated,
               OLCompose(@(1).val, @(2).val._children[1]))]),
    Reduction_ScatHUnion := ARule(OLCompose, [ @(1, Reduction), @(2, ScatHUnion, e->e.n=1) ], e->[]),
    ISumReduction_PointWise := ARule(OLCompose, [ @(1, ISumReduction), @(2, PointWise) ],
        e -> [ ISumReduction(@(1).val.var, @(1).val.domain, @(1).val.op, @(1).val.idval, @(1).val.isSaturated,
                OLCompose(@(1).val._children[1], @(2).val))]),
    eT_Pointwise := ARule(OLCompose, [@(1, eT, e->IsVar(e.base)), @(2, PointWise)],
        e -> let(i := Ind(1),
            [PointWise(1, Lambda(
                SubstVars(Copy(@(2).val.op.vars), rec((@(2).val.op.vars[2].id):=i)),
                SubstVars(Copy(@(2).val.op.expr), rec((@(2).val.op.vars[2].id):=@(1).val.base)))), @(1).val])),
    ISumXXX_YYY := ARule(OLCompose, [ @(1, [ISumUnion, ISumReduction]), @(2, [GathH, PointWise, Induction, eT]) ],
        e -> [ CopyFields(@(1).val, rec(
            _children :=  List(@(1).val._children, c -> OLCompose(c, @(2).val)),
            dimensions := [Rows(@(1).val), Cols(@(2).val)] )) ]),
    GathH_GathH := ARule(OLCompose, [@(1, GathH), @(2, GathH)],
        e -> [GathH(@(2).val.N, @(1).val.n, @(1).val.base+@(2).val.base, @(1).val.stride*@(2).val.stride)]),
    PointWise_BinOp := ARule(OLCompose, [@(1, PointWise, e->e.N=1), @(2, BinOp, e->e.N=1)],
        e -> [BinOp(1,Lambda(@(2).val.op.vars, @(1).val.op.at(@(2).val.op.expr, V(0))))]),
    eT_Induction := ARule(OLCompose, [@(1, eT), @(2, Induction)],
        e -> [Inductor(@(2).val.N, @(1).val.base, @(2).val.op, @(2).val.initval)]),
    ScatHUnion_ScatHUnion := ARule(OLCompose, [@(1, ScatHUnion), @(2, ScatHUnion)],
        e -> [ScatHUnion(@1.val.N, @(2).val.n, @(1).val.base+@(2).val.base, @(2).val.stride)] )
));

Class(ToVectors, HierarchicalVisitor, rec(
    __call__ := meth(arg)
    local res;
        res := ApplyFunc(arg[1].visit, arg{[2..Length(arg)]});
       return res;
    end,
    PointWise := (self, o, opts) >> let(
      new_i := Ind(4), new_var := var.fresh_t("r", TReal), times := (o.N - Mod(o.N, 4))/4,
      DirectSum(
        List([0..times-1], i->PointWise(4, Lambda([new_var, new_i], new_var*i)) ),
	PointWise(Mod(o.N, 4), o.op)
      )
    )
));

max_value := var("max_val", TReal);
max_err := var("max_err", TReal);
Class(RulesErrorHCOL, RuleSet, rec(inType := "iCode", outType := "iCode"));
RewriteRules(RulesErrorHCOL, rec(
  assign_const := Rule([@(1, assign), @(2, var), @(3, Value, v->v<>0)],
					e->assign(@(2).val, @(3).val+@(3).val*max_err) ),
  assign_nth := Rule([@(1, assign), @(2, var), @(3, nth)],
					e->assign(@(2).val, max_value+max_value*max_err) )
));

Class(RulesCodeHCOL, RuleSet, rec(inType := "iCode", outType := "iCode"));
RewriteRules(RulesCodeHCOL, rec(
#only assignments to state but state is not used

	simplify_state := Rule(@@(1, decl, (e, cx)->(
	  let(
	  Length(Collect(e.cmd, decl))=0 and
	  IsBound(cx.opts.state) and
	  Length(Collect(e.cmd, @(1,nth).cond(j->IsVar(j.loc) and j.loc.id=cx.opts.state.id))) > 0 and
	  Length(Collect(e.cmd, @(1,nth).cond(j->IsVar(j.loc) and j.loc.id=cx.opts.state.id)))=Length(Collect(e.cmd, @(1, assign).cond(i->ObjId(i.loc)=nth and i.loc.loc.id = cx.opts.state.id))) ))),
	  (e, cx)-> let(Print("Remove State:\n", @@(1).val, "\n"), SubstTopDown(@@(1).val, [@(2, assign), @(3, nth).cond(i->i.loc.id=cx.opts.state.id), @(4)], j->skip()))),

    chain_chain := Rule(@(1, chain, e->ObjId(e.cmds[1])=chain),
        e -> let(chain(Concat(@(1).val.cmds[1].cmds, Drop(@(1).val.cmds, 1))))),

	#replaces TArray of length one with a scalar variable
    scalarize1 := Rule(@(1, decl, e->Length(Filtered(e.vars, i->ObjId(i.t)=TArray and i.t.size=1))>=1),
        e->let(
			oldvar := Filtered(e.vars, i->ObjId(i.t)=TArray and i.t.size=1)[1],
			Print("Scalarizing ", oldvar, "\n"),

			tmp := Collect(@(1).val, [@(5, assign), @(6).cond(ee->IsNth(ee) and ee.loc.id=oldvar.id), @(7)]),
			newType := UnifyTypes(List(tmp, i->i.exp.t)),
			newvar := var.fresh_t("s", newType),
			Print("Scalarizing ", oldvar, " with ", newvar, " of type: ", newType, "\n"),

            decl(
                SubstVars(@(1).val.vars, rec((oldvar.id):=newvar)),
                SubstTopDown(@(1).val.cmd, [@(2, nth), @(3).cond(e->IsVar(e) and e.id=oldvar.id), @(4).cond(e->IsValue(e) and e.v = 0)],
                    i->newvar)
            ))),

	#replace each fixed index array reference (i.e. nth(X, i), where i is a constant value) with a variable
	#This can only be done if there are no array references to X that is variable (this is to prevent potential aliasing)
    scalarize_const := Rule(@(1, decl, e->Length(Filtered(Filtered(e.vars, j->IsArray(j.t)), k->ForAll(Collect(e.cmd, [@(2, nth), @(3, var, u->u.id=k.id), @(4)]), j->IsValue(j.idx))))>=1),
        e->let(
			oldvar := Filtered(Filtered(e.vars, j->IsArray(j.t)),   #find all array variables
	           k->ForAll(Collect(e.cmd, [@(2, nth), @(3, var, u->u.id=k.id), @(4)]), j->IsValue(j.idx)))[1],   #make sure only const index

            newvars := List([1..oldvar.t.size], j->var.fresh_t("q", oldvar.t.t)),
            cmmd := @(1).val.cmd,
            lst := List([0..Length(newvars)-1], i->[i, newvars[i+1]]),
            f := (c,l)->SubstTopDown(c, [@(2, nth), @(3).cond(k->IsVar(k) and k.id=oldvar.id), @(4).cond(j->IsValue(j) and j.v = l[1])], i->l[2]),
            Print("Scalarize_const: ",oldvar,":: ", newvars, "\n"),
            decl(Concat(Filtered(@(1).val.vars, j->j.id <> oldvar.id), newvars), FoldL(lst, f, cmmd))
        )),

    chain_xyz_chain := ARule(chain, [@(1), @(2, chain)],
        e->[chain(Concat([@(1).val], @(2).val.cmds))]),
    chain_chain_xyz := ARule(chain, [@(1, chain), @(2)],
        e->[chain(Concat(@(1).val.cmds), [@(2).val])]),

    chain_creturn := ARule(chain, [@(1, creturn), @(2)],
        e->[chain(Concat([@(2).val], [@(1).val]))] ),

    loop_pull_cond_1 := Rule([@(1, loop), @(2), @(3), [@(4, chain),  [@(5, assign), @(6,var),
            [@(7,cond), @(8,eq, e->e.args[1]=@(2).val and e.args[2].v=0), @(9), @(10)]],...]],
        e->let(
	   chain(assign(@(6).val, @(9).val), loop(@(1).val.var, @(1).val.range, chain([assign(@(6).val, @(10).val)]::Drop(@(4).val.cmds,1)))))),

	   copyprop_var_val := Rule(@(1, decl, e->let(
		   #only propagate values if the variable is SSA and only if it is within a basic block
		   cnds := Collect(e.cmd, [assign, @(2, var, ee->ee in e.vars), @(3, Value)]),
		   actual_cnds := Filtered(cnds, i->Length(Collect(e.cmd, [assign, @(8, var, ee->ee.id = i.loc.id), @(9)]))=1),
		   #Print("CopyProp: ", cnds, "\n\n", e.cmd, "\n\n", actual_cnds, "\n"),
		   Length(actual_cnds) >= 1 and Length(Collect(e.cmd, [assign, @(4, var, ee->ee.id = actual_cnds[1].loc.id), @(5)]))=1)),
		e->let(
			all_assigned_value := Collect(@(1).val.cmd, [assign, @(6, var, ee->ee in e.vars), @(7, Value)]),
			assigned_var_value := Filtered(all_assigned_value, i->Length(Collect(@(1).val.cmd, [assign, @(8, var, ee->ee.id = i.loc.id), @(9)]))=1)[1],
		   #Print("CopyProp_value: ", assigned_var_value, "\n"),
			decl(Filtered(@(1).val.vars, j->j.id<>assigned_var_value.loc.id),
				 SubstVars(@(1).val.cmd, rec((assigned_var_value.loc.id) := assigned_var_value.exp)))
		)
       ),

	   copyprop_var_var := Rule(@(1, decl, e->let(
		   #only propagate var if the assigned variable is SSA
           cnds := Collect(e.cmd, [assign, @(2, var, ee->ee in e.vars), @(3, var, e->e.id<>@(2).val.id)]),
           Length(cnds)>=1 and cnds[1].loc in e.vars and Length(Collect(e.cmd, [assign, @(4, var, ee->ee.id=cnds[1].loc.id), @(5)]))=1)),
        e->let(
			   acnds := Collect(@(1).val.cmd, [assign, @(7, var, ee->ee in e.vars), @(8, var, e->e.id<>@(7).val.id)]),
			   atests := Collect(e.cmd, [assign, @(9, var, ee->ee.id=acnds[1].loc.id), @(5)]),
			   a := atests[1],
			   decl(Filtered(@(1).val.vars, j->j <>a.loc), SubstVars(@(1).val.cmd, rec((a.loc.id) := a.exp))))),

    drop_selfassign_var := ARule(chain, [@(1, assign, e->IsVar(e.loc) and IsVar(e.exp) and e.exp=e.loc), @(2)],
        e->let(
		[@(2).val])),

    drop_selfassign2 := ARule(chain, [@(1, assign, e->e.exp=e.loc), @(2)],
        e->[@(2).val]),

	#this drops unused variables i.e. variables that are declared but not used.
    drop_unused_var := Rule(@(1, decl, e-> let(v := Collect(e.cmd, var), Filtered(e.vars, k->not k in v))<> []),
        e->let(v := Collect(@(1).val.cmd, var),
            decl(Filtered(@(1).val.vars, k->k in v), @(1).val.cmd))),

	drop_unused_assign := Rule(@(1, decl, e->
			Filtered(List(e.vars, v->Collect(e.cmd, v)), l->Length(l)=1)<>[]) ,
	e->let(
      unused := Filtered(List(e.vars, v->Collect(e.cmd, v)), l->Length(l)=1)[1],
	  empty_list := List([1..Length(unused)], i->skip()),
	  target_list := Filtered(Collect(@(1).val.cmd, assign), a->a.loc in unused),
	  Print("Unused assignements:\n", target_list, "\n\n"),
	  decl(Filtered(@(1).val.vars, k->not k in unused), SubstList(@(1).val.cmd, target_list, empty_list))
	)),

    loop_decl := Rule(@(1, loop, e->ObjId(e.cmd)=decl and Length(Collect(e.cmd, decl))=1 and Length(Collect(e.cmd, chain))<=1),
        e -> decl(@(1).val.cmd.vars, loop(@(1).val.var, @(1).val.range, @(1).val.cmd.cmd))),

    chain_xyz_decl := ARule(chain, [@(1), @(2, decl)],
        e->[decl(@(2).val.vars, chain(@(1).val, @(2).val.cmd))]),

    chain_decl := Rule(@(1, chain, e->ObjId(e.cmds[1])=decl),
        e -> decl(@(1).val.cmds[1].vars, chain(Concat([@(1).val.cmds[1].cmd], Drop(@(1).val.cmds, 1))))),
    decl_decl := Rule(@(1, decl, e->ObjId(e.cmd)=decl),
        e -> decl(Set(Concat(@(1).val.vars, @(1).val.cmd.vars)), @(1).val.cmd.cmd)),
));

Class(RulesHCOLnoAbsMax, RuleSet, rec(inType := "iCode", outType := "iCode"));
RewriteRules(RulesHCOLnoAbsMax, rec(
    abs_cond := Rule([@(1, assign), @(2), @(3, abs)],
        e->let(s := var.fresh_t("w", @(3).val.t), decl(s, chain(assign(s, @(3).val.args[1]), assign(@(1).val.loc, cond(geq(s, V(0)), s, neg(s))))))),
    max_cond := Rule(@(1, max),
        e->cond(geq(@(1).val.args[1], @(1).val.args[2]), @(1).val.args[1], @(1).val.args[2])),
));

Class(RulesUnrollHCOL, RuleSet, RulesStrengthReduce, RulesCodeHCOL, rec(inType := "iCode", outType := "iCode"));
RewriteRules(RulesUnrollHCOL, rec(
    unroll_wo_decl := Rule(@@(1, [loopn, loop],
       (e, cx) -> Length(Collect(e.cmd, loop))=0 and
                  Length(Collect(e.cmd, decl))=0 and
	          e.var.range <= When(IsBound(cx.opts) and IsBound(cx.opts.globalUnrolling),
		       	                                           cx.opts.globalUnrolling,
                                              5)),
       (e,cx)->let(i := @@(1).val.var,
	               rng := @@(1).val.range,
		       vvars := Filtered(Collect(@@(1).val.cmd, var), j->j.t in [TReal, TDouble]),

           chain(
	      List(rng, j->let(
		    ssvars := FoldL([[i, V(j)]], (b, a)->CopyFields(rec((a[1].id):= a[2]), b), rec()),
		    new_cmd := SubstVars(Copy(@@(1).val.cmd), ssvars),
 		    new_cmd ))
    ))) ,
    unroll_w_decl := Rule(@@(1, [loop, loopn],
           (e,cx)->Length(Collect(e.cmd, loop))=0 and
				   Length(Collect(e.cmd, decl))=1 and
	           e.var.range <= When(IsBound(cx.opts) and IsBound(cx.opts.globalUnrolling),
           	                                          cx.opts.globalUnrolling,
                                              5)),
           (e,cx)->let(i := @@(1).val.var,
	               rng := @@(1).val.range,
		       vvars := Filtered(Collect(@@(1).val.cmd, var), j->j.t in [TReal, TDouble]),

		       localvar := Collect(@@(1).val.cmd, decl)[1].vars,

           chain(
	      List(rng, j->let(

	      		        nvars :=FoldL(localvar, (b,a)->Concat([[a, var.fresh_t("u", a.t)]], b), [[i, V(j)]] ),
	      		        ssvars := FoldL(nvars, (b, a)->CopyFields(rec((a[1].id):= a[2]) , b), rec()),
						svars := FoldL(nvars, (b, a)->Concat( When(IsVar(a[2]) or a[2] in b, [a[2]], []), b), []),

				new_cmd := SubstVars(Copy(@@(1).val.cmd.cmd), ssvars),

				decl(svars, new_cmd) )

                                #decl(svars, new_cmd ) )
             ))
	  )
   )
));

RewriteRules(RulesStrengthReduce, rec(
	addsub00 := Rule([@(1, addsub_2x64f), @(2).cond(e->Cond(IsValue(e), isValueZero(e), e=0)), _0],
		e->@(2).val)
));

Class(RulesUnifyType, RuleSet, rec(inType := "iCode", outType := "iCode"));
RewriteRules(RulesUnifyType, rec(
    unify_value := Rule(@(1,assign, e-> IsValue(e.exp) and IsVar(e.loc) and e.loc.t <> e.exp.t),
                        e->let(
                            newtype := UnifyTypes([@(1).val.exp.t, e.loc.t]),
                            Print("Unifying Type from a value to ", newtype,"\n"),
                            #if the newtype is the same as the loc, then need to cast the value
                            #else error, upcast the type of the variable since we are losing info.
                            Cond(newtype=e.loc.t,
                            assign(@(1).val.loc, tcast(newtype, @(1).val.exp)),
                            let(
                            update_type(@(1).val.loc, newtype),
                            assign(@(1).val.loc, @(1).val.exp)
                            ))
                        )),
	unify_vartype := Rule(@(1,assign,e->IsVar(e.loc) and not IsValue(e.exp) and e.loc.t <> e.exp.computeType()),
	   e->let(
		  oldtype := @(1).val.exp.t,
		  newtype := @(1).val.exp.computeType(),
		  Print("Unifying Type from ", oldtype, " to ", newtype,"\n"),
		  update_type(@(1).val.loc, newtype),
		  assign(@(1).val.loc, @(1).val.exp)
	   )),
));
RulesTypeHCOL := CopyFields(MergedRuleSet(RulesUnifyType),
    rec(inType:="iCode", outType := "iCode"));

RulesOrigCodeUnrollHCOL := CopyFields(MergedRuleSet(RulesUnrollHCOL, RulesStrengthReduce, RulesCodeHCOL, RulesHCOLnoAbsMax),
    rec(inType:="iCode", outType := "iCode"));

RulesCodeUnrollHCOL := CopyFields(MergedRuleSet( RulesUnrollHCOL, RulesStrengthReduce,  RulesCodeHCOL, RulesHCOLnoAbsMax),
    rec(inType:="iCode", outType := "iCode"));

RulesCodeNoUnrollHCOL := CopyFields(MergedRuleSet(RulesStrengthReduce, RulesCodeHCOL, RulesHCOLnoAbsMax),
    rec(inType:="iCode", outType := "iCode"));

#changes behaviour so that Spiral outputs abs(a) and max(a,b) C-code instead of "(a >= b) ? a : b"
RulesCodeUnrollHCOLuseAbsMaxSR := CopyFields(MergedRuleSet(RulesUnrollHCOL, RulesStrengthReduce, RulesCodeHCOL),
    rec(inType:="iCode", outType := "iCode"));

max_stars_repo_name: spiral-software/spiral-package-hcol

max_issues_repo_name: spiral-software/spiral-package-hcol

max_forks_repo_name: spiral-software/spiral-package-hcol


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

Declare(DFT_Rader, DFT_GoodThomas);

_accurate_dft := vec -> List(TransposedMat(MatSPL(transforms.DFT(Length(vec)))*TransposedMat([vec]))[1], ComplexCyc);

Class(fdft, FuncExp, rec(
    computeType := self >> self.args[1].t
));

Class(fDFT, DiagFunc, rec(
    abbrevs := [ func -> Checked(IsFunction(func), [func]) ],
    range  := self >> TComplex,
    domain := self >> self.params[1].domain(),

    lambda := self >> fdft(self.params[1].lambda()),

    tolist := self >> let(
	elts := List(self.params[1].tolist(), x->Complex(EvalScalar(x))),
	ComplexFFT(elts))
));

NewRulesFor(DFT, rec(
    #F DFT_Base: DFT_2 = F_2
    #F
    DFT_Base := rec(
        info             := "DFT_2 -> F_2",
        forTransposition := false,
        applicable       := nt -> nt.params[1] = 2 and not nt.hasTags(),
        apply            := (nt, C, cnt) -> F(2)
    ),

    DFT_Base1 := rec(
        info             := "DFT_1 -> I(1)",
        forTransposition := false,
        applicable       := nt -> nt.params[1] = 1 and not nt.hasTags(),
        apply            := (nt, C, cnt) -> I(1)
    ),

    #F DFT_Canonize: DFT_n_k = perm * DFT_n
    #F
    DFT_Canonize := rec(
       switch           := false,
       info             := "DFT_n_k -> perm * DFT_n",
       forTransposition := false,
       applicable       := nt -> nt.params[1] > 2 and nt.params[2] <> 1 and not nt.hasTags(),
       children         := nt -> [[ DFT(nt.params[1], 1) ]],
       apply            := (nt, C, cnt) -> OddStride(nt.params[1], nt.params[2]) * C[1]
    ),

    #F DFT_CT: 1965
    #F   General Cooley-Tukey Rule
    #F   DFT_n = (DFT_n/d tensor I_d) * diag * (I_n/d tensor F_d) * perm
    #F
    #F Cooley/Tukey:
    #F   An Algorithm for the Machine Calculation of Complex Fourier Series.
    #F   Mathematics of Computation, Vol. 19, 1965, pp. 297--301.
    #F
    DFT_CT := rec(
        info          := "DFT(mn,k) -> DFT(m, k%m), DFT(n, k%n)",

        maxSize       := false,
        forcePrimeFactor := false,

        applicable := (self, nt) >> nt.params[1] > 2
            and not nt.hasTags()
            and (self.maxSize=false or nt.params[1] <= self.maxSize)
            and not IsPrime(nt.params[1])
            and When(self.forcePrimeFactor, not DFT_GoodThomas.applicable(nt), true),

        children  := nt -> Map2(DivisorPairs(nt.params[1]),
            (m,n) -> [ DFT(m, nt.params[2] mod m), DFT(n, nt.params[2] mod n) ]
        ),

        apply := (nt, C, cnt) -> let(mn := nt.params[1], m := Rows(C[1]), n := Rows(C[2]),
            Tensor(C[1], I(n)) *
            Diag(fPrecompute(Tw1(mn, n, nt.params[2]))) *
            Tensor(I(m), C[2]) *
            L(mn, m)
        )
    ),

    DFT_CT_Mincost := rec(
        maxSize := false,
        applicable := (self, nt) >> let(N := nt.params[1],
	    N > 2 and not nt.hasTags() and IsEvenInt(N) and
            (self.maxSize = false or nt.params[1] <= self.maxSize)),

        children  := nt -> let(N := nt.params[1], k := nt.params[2], Map2(
	    Filtered(DivisorPairs(N), divpair -> IsEvenInt(divpair[2])),
            (m,n) -> [ DFT(m, k mod m), DFT3(m, k mod (2*m)), DFT(n, k mod n) ])),

        apply := (nt, C, cnt) -> let(
	    N := nt.params[1], m := Rows(C[1]), n := Rows(C[3]), k := nt.params[2],
            hf := (n-2)/2,     i := Ind(hf),    j := Ind(hf),

            SUM(        _hhi(N,  m, 0,      n,  C[1]), # first iteration has no twiddles
              When(hf=0, [],
		ISum(i, _hhi(N,  m, i+1,    n,  C[1] * Diag(fPrecompute(Twid(N,m,k,0,0,i+1)))))),
                        _hhi(N,  m, hf+1,   n,  C[2]), # middle is a DFT3
	      When(hf=0, [],
		ISum(j, _hhi(N,  m, j+hf+2, n,  C[1] * Diag(fPrecompute(Twid(N,m,k,0,0,j+2+hf))))))
            ) *
            Tensor(I(m), C[3]) *
            Tr(n, m)
        )
    ),

    #F DFT_CosSinSplit:
    #F
    #F   DFT_n = CosDFT_n + i * SinDFT_n
    #F
    #F This rule has been restricted in applicability to small size
    #F (see config.g), Connects to Vetterli rules for CosDFT/SinDFT.
    #F
    DFT_CosSinSplit := rec(
        info             := "DFT_n -> CosDFT_n, SinDFT_n",
        forTransposition := false,
        maxSize := 64,

        applicable := (self, nt) >> let(N:=nt.params[1], rot:=nt.params[2],
            N > 2 and N <= self.maxSize and Is2Power(N) and AbsInt(rot)=1
	    and not nt.hasTags()),

        children := nt -> [[ CosDFT(nt.params[1]), SinDFT(nt.params[1]) ]],

        apply := (nt, C, cnt) -> let(N := nt.params[1], rot := nt.params[2],
	    Tensor(Mat([[1,1]]), I(N)) *
            VStack(C[1], (E(4)^rot) * C[2])
	)
    ),

    #F DFT_Rader: Rader's Algorithm for Prime size DFT
    #F
    #F   DFT(p) -> P(g)' * (1 dirsum DFT(p-1)) * Tp * (1 dirsum DFT(p-1)) * P(g)",
    #F   P(g) = perm
    #F   Tp   = [[1,1],[1,-1/(p-1)]] dirsum diag
    #F
    DFT_Rader := rec(
        info             := "DFT(p) -> P(g)' * DFT(p-1) * Tp ",
        forTransposition := false,
        minSize          := 3,
        maxSize          := 64,
        accurate         := true,

        raderDiag := (self, N, k, root) >> let(
	    dft := When(self.accurate, _accurate_dft, vec->ComplexFFT(List(vec,ComplexCyc))),
            SubList(dft(List([0..N-2], x -> E(N)^(k*root^x mod N))), 2) / (N-1)),

        raderMid := (self, N, k, root) >>
                     DirectSum(Mat([[1, 1], [1, -1/(N-1)]]),
                       Diag(FData(self.raderDiag(N, k, root)))),

        applicable := (self, nt) >> let(n:=nt.params[1],
            n >= self.minSize and n <= self.maxSize and IsPrime(EvalScalar(n))
	    and not nt.hasTags()),

        children      := nt -> [[ DFT(nt.params[1]-1, -1), DFT(nt.params[1]-1, -1) ]],

        apply := (self, nt, C, cnt) >> let(
	    N := EvalScalar(nt.params[1]), k := EvalScalar(nt.params[2]), root := PrimitiveRootMod(N),
            RR(N, 1, root).transpose() *
            DirectSum(I(1), C[1]) *
            self.raderMid(N,k,root) *
            DirectSum(I(1), C[2]) *
            RR(N, 1, root)
        )
    ),

    DFT_RealRader := CopyFields(~.DFT_Rader, rec(
        children      := nt -> [[ PRDFT1(nt.params[1]-1).transpose(), PRDFT1(nt.params[1]-1) ]],

        realRaderDiag := (self, N, k, root) >> let(cdiag := self.raderDiag(N, k, root),
            # cdiag is of length N-2.
            # cdiag = sublist( DFT_{N-1} * omega_vector, [2..N-1])
            # now we want to obtain
            #   rdiag = sublist( RDFT_N * omega_vector, [2..N-1]), we can do this using
            #   the property RDFT = (1 dirsum (I tensor [[1, 1], [-j, j]]) dirsum 1) * LIJ * DFT
            #   which means that rdiag = (1 dirsum (I tensor [[1, 1], [-j, j]]) dirsum 1) * LIJ * cdiag
            #   the code below does this ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ plus pads a zero (for PRDFT)
            #   and also the following trick.

            #  The expected rdiag steady state is [[a, b], [b, -a]], we convert this to rotations (=RCDiag)
            #  by applying J(2), ie we get [[b, -a], [a, b]], which would be a rotation by b + E(4)*a,
            #  but in this case both a and b are complex, so its not. We still use RCDiag though.
            #
            Concatenation(
                ConcatList([1..(N-3)/2], i -> Reversed([cdiag[i]+cdiag[N-1-i], -E(4)*(cdiag[i]-cdiag[N-1-i])])), # DFT->RDFT
                [0, cdiag[(N-3)/2+1]])),

        raderMid := (self, N, k, root) >>
            DirectSum(Mat([[1, 1, 0], [1, -1/(N-1), 0], [0, 0, 0]]), # zero padding for PRDFT
                Tensor(I((N-1)/2), J(2))*
                RCDiag(FData(self.realRaderDiag(N, k, root))).toloop()),
                # RCDiag with complex entries, max weird :) !
                # .toloop() is used because in vector code RCDiag is
                # converted to VRCDiag which is confused by complex entries AFAIK
        apply := (self,nt,C,cnt) >> let(
        N := EvalScalar(nt.params[1]), k := EvalScalar(nt.params[2]), root := PrimitiveRootMod(N),
            RR(N, 1, root).transpose() *
            DirectSum(I(1), C[1]) *
            self.raderMid(N,k,root) *
            DirectSum(I(1), C[2]) *
            RR(N, 1, root))
    )),

    #F DFT_Bluestein : Convert FFT to Toeplitz matrix, then embed in a larger Circulant
    #F
    #F    DFT_n -> diag * Toeplitz * diag
    #F
    # DFT on inputs, some of which are 0:
    #
    DFT_Bluestein := rec(
        info             := "DFT_n -> diag * Toeplitz * diag",
        forTransposition := false,
        minSize := 3,
        maxSize := false,
        switch := false,

        # Applicable only for non-powers of 2 to prevent cycles
        applicable := (self, nt) >> let(N := nt.params[1],
            IsSymbolic(N) or
	    (N >= self.minSize and When(IsInt(self.maxSize), N <= self.maxSize, true) and not Is2Power(N)
		and not nt.hasTags())),

	# NOTE: extend this to be a true degree of freedom
	freedoms := t -> let(min := 2*t.params[1]+1,
	    [ integersBetween(min, 4 * pow(2, floor(log(t.params[1])/log(2)))) ]),

        child := (t, fr) -> let(
	    N := t.params[1], Ncirc := fr[1],
            [ DFT(Ncirc, -1).withTags(t.getTags()), DFT(Ncirc, 1).withTags(t.getTags()) ]),

        diag := (N, k) -> let(i:=Ind(N), Lambda(i, omegapi(fdiv(i^2 * k, N)))),

	circ := (N, k, circsize) -> let(
	    j1 := Ind(N),
	    j2 := Ind(N-1),
	    diagDirsum(
		Lambda(j1, fdiv(omegapi(fdiv(-j1^2*k, N)), circsize)),
		fConst(TReal, circsize-2*N+1, 0),
		Lambda(j2, fdiv(omegapi(fdiv(-(N-j2-1)^2*k, N)), circsize)))),

        apply := (self, nt, C, cnt) >> let(
            N := nt.params[1],
            k := nt.params[2],
            Ncirc := Rows(C[1]),
            diag := fPrecompute(self.diag(N,k)),
	    pad  := spiral.paradigms.common.TRDiag(Ncirc, N, diag).withTags(nt.getTags()),

	    #Diag(diag) * Gath(fCompose(fAdd(Ncirc, N, 0))).toloop(1) *
	    pad.transpose() *
	    Inplace(Grp(C[1] * Diag(fPrecompute(fDFT(self.circ(N, k, Ncirc)))))) *
	    Inplace(C[2]) *
	    pad
	    #Scat(fAdd(Ncirc, N, 0)).toloop(1) * Diag(diag) ##, O(Ncirc-N, N))
	)
    ),

    #F DFT_GoodThomas : Prime Factor FFT
    #F
    #F     DFT_n*k -> perm * (DFT_n_a tensor DFT_k_b) * perm
    #F     when gcd(n,k) = 1
    #F
    DFT_GoodThomas := rec(
        info             := "DFT_n*k -> perm * (DFT_n_a tensor DFT_k_b) * perm",
        forTransposition := false,
        applicable     := (self, nt) >> let(N := nt.params[1],
	    N > 2 and DivisorPairsRP(N) <> [] and N <= self.maxSize and
	    not nt.hasTags()),

        maxSize := 64,

        children := nt -> let(mn := nt.params[1], k:= nt.params[2], Map2(DivisorPairsRP(mn),
            (m,n) -> [ DFT(m, k*n mod m), DFT(n, k*m mod n) ])),

        apply := (nt, C, cnt) -> let(
            r  := Rows(C[1]),
            s  := Rows(C[2]),
            alpha := 1 / s mod r,
            beta  := 1 / r mod s,
            CRT(r,s,1,1).transpose() * Tensor(C[1], C[2]) * CRT(r,s,1,1)
        )
    ),

    #F DFT_PFA_SUMS : Prime Factor FFT using Sigma-SPL and cyclic shifts
    #F
    #F     DFT_n*k -> perm * (DFT_n_a tensor DFT_k_b) * perm
    #F     when gcd(n,k) = 1
    #F
    DFT_PFA_SUMS := CopyFields(~.DFT_GoodThomas, rec(
        apply := (nt, C, cnt) -> let(
	    r := Rows(C[1]), s := Rows(C[2]),
	    j := Ind(s), k := Ind(r),
	    alpha := (1/s) mod r, beta := (1/r) mod s,
	    ISum(j,
                Scat(fTensor(fId(r), fBase(j))) *
                (C[1] ^ Z(r, -alpha*j)) *
                Gath(fTensor(fId(r), fBase(j)))) *
	    ISum(k,
                Scat(fTensor(fId(s), fBase(k))) *
                (C[2] ^ Z(s, -beta*k)) *
                Gath(fTensor(fId(s), fBase(k))))
	)
    )),

    DFT_PFA_RaderComb := rec(
        forTransposition := false,
        applicable := nt -> let(N := nt.params[1], divs := DivisorPairsRP(N),
            N > 2 and Length(divs)=2 and N mod 2 = 1 and
            IsPrime(divs[1][1]) and IsPrime(divs[1][2])
	    and not nt.hasTags()),

        children := nt -> let(
	    N:=nt.params[1], k:=nt.params[2], divs := DivisorPairsRP(nt.params[1]),
             r := divs[1][1], s := divs[1][2],
                 [ [DFT(r-1,-1), DFT(r-1,-1), DFT(s-1,-1), DFT(s-1,-1)],
		   [DFT(s-1,-1), DFT(s-1,-1), DFT(r-1,-1), DFT(r-1,-1)] ]),

        apply := (nt, C, cnt) -> let(
            r := Rows(C[1])+1,
            s := Rows(C[3])+1,
            er := ChineseRem([r,s], [1,0]),
            es := ChineseRem([r,s], [0,1]),
            dr := DFT_Rader.raderMid(r, nt.params[2]*s), # * er/s mod r),
            ds := DFT_Rader.raderMid(s, nt.params[2]*r), # * es/r mod s),

            CRT(r,s,1,1).transpose() *
            Tensor(RR(r).transpose(), RR(s).transpose()) *
            Tensor(
                DirectSum(I(1), C[1]) * dr * DirectSum(I(1), C[2]),
                DirectSum(I(1), C[3]) * ds * DirectSum(I(1), C[4])
            ) *
            Tensor(RR(r), RR(s)) *
            CRT(r,s,1,1)
        )
    ),

    #F DFT_SplitRadix: 1984
    #F
    #F DFT_n = B * (DFT_n/2 dirsum DFT_n/4 dirsum DFT_n/4) * perm
    #F
    #F B = (DFT_2 tensor I_n/2) * S * diag
    #F S = (I_n/2 dirsum (diag([1,E(4)] tensor I_n/4))
    #F
    #F Duhamel, Pierre:
    #F   Split Radix FFT Algorithm, Electronics Letters, Vol. 20, No. 1,
    #F   pp. 14--16, 1984
    #F
    DFT_SplitRadix := rec(
        info             := "DFT_n -> DFT_n/2, DFT_n/4, DFT_n/4",
        forTransposition := true,
        maxSize := 64,

        applicable := (self, nt) >> let(N := nt.params[1],
            N >= 8 and N mod 4 = 0 and N <= self.maxSize and not nt.hasTags()
        ),

        children := nt -> let(N := nt.params[1], w := nt.params[2],
            [[ DFT(N/2, w), DFT(N/4, w) ]]
        ),

        apply := (nt, C, cnt) -> let(N := nt.params[1], w := nt.params[2],
            Tensor(F(2), I(N/2)) *
            DirectSum(
                C[1],
                Tensor(Diag([1, E(4)^w]) * F(2), I(N/4)) *
                Diag(fPrecompute(fCompose(Tw3(N/2, 2, w), L(N/2, 2)))) *
                Tensor(I(2), C[2])*L(N/2,2)
            ) *
            L(N, 2)
        )
    ),

    #F DFT_DCT1andDST1: 1984
    #F
    #F   DFT_n = blocks * (DCT1_(n/2+1) dirsum DST1_(n/2-1)) * blocks
    #F
    #F   Wang:
    #F     Fast Algorithms for the Discrete W Transform and the
    #F     Discrete Fourier Transform.
    #F     IEEE Trans. on ASSP, 1984, pp. 803--814
    #F   Britanak/Rao:
    #F     The fast generalized discrete Fourier transforms: A unified
    #F     approach to the discrete sinusoidal transforms computation.
    #F     Signal Processing 79, 1999, pp. 135--150
    #F
    DFT_DCT1andDST1 := rec(
        info    := "DFT_n -> DCT1_(n/2+1), DST1_(n/2-1)",
        maxSize := 64,

        applicable := (self, nt) >> let(N:=nt.params[1], k:=nt.params[2],
            N > 2 and N <= self.maxSize and Is2Power(N) and AbsInt(k)=1) and not nt.hasTags(),

        children := nt ->
	    [[ DCT1( nt.params[1]/2 + 1 ), DST1( nt.params[1]/2 - 1 ) ]],

        apply := (nt, C, cnt) -> let(N:=nt.params[1], w:=nt.params[2],
            DirectSum(I(1), blocks4(N - 1, w)) *
            DirectSum(C[1], C[2] ^ J(N/2 - 1)) *
            DirectSum(I(1), blocks1(N - 1))
        )
    )
));