text:
# A Taggable operator (see tags.gi) is translated into an analogous SIMT operator that carries dim tag info

Class(SIMTTensor, Tensor, rec(

    abbrevs := [ (simt_dim, spl) -> When(IsASIMTDim(simt_dim) and ObjId(spl) = Tensor, [Flat([simt_dim, spl._children])], [[simt_dim, spl]] )  ]
                :: Tensor.abbrevs,

    new := meth(self, arg)
        local simt_dim, L;
        Constraint(IsList(arg) and Length(arg) >= 2);
        simt_dim := arg[1];
        L := arg{[2..Length(arg)]};
#        L:=Filtered(L,x->x<>I(1));
        if (Length(L)=0) then return I(1); fi;
        if Length(L)=1 then return L[1]; fi;

        return SPL(WithBases( self,
                        rec( simt_dim := simt_dim,
                             _children := L,
                             dimensions := [ Product(L, t -> t.dims()[1]),
                                             Product(L, t -> t.dims()[2]) ] )));
    end,

    from_rChildren := (self, rch) >> ApplyFunc(ObjId(self), [self.simt_dim]::rch).appendAobj(self)

    ));

Class(SIMTISum, ISum, rec(

    abbrevs := [ (simt_dim, spl) -> When(IsASIMTDim(simt_dim) and ObjId(spl) = ISum, [simt_dim, spl.var, spl.domain, spl._children[1]], [simt_dim, spl] )  ],

    new := meth(self, simt_dim, var, domain, expr)
        local res;
        Constraint(IsSPL(expr));
        # if domain is an integer (not symbolic) it must be positive
        Constraint(not IsInt(domain) or domain >= 0);
        var.isLoopIndex := true;
        res := SPL(WithBases(self, rec(_children := [expr], simt_dim := simt_dim, var := var, domain := domain)));
        res.dimensions := res.dims();
        return res;
    end,

    print := (self, i, is) >> Print(
        self.name, "(", self.simt_dim, ", ", self.var, ", ", self.domain, ",\n",
        Blanks(i+is), self._children[1].print(i+is, is), "\n",
        Blanks(i), ")", self.printA(),
        When(IsBound(self._setDims), Print(".overrideDims(", self._setDims, ")"), Print(""))
    )

    ));

Class(SIMTSUM, SUM, rec(

    new := meth(self, arg)
        local dims, simt_dim, L;
        Constraint(IsList(arg) and Length(arg) >= 2);
        simt_dim := arg[1];
        L := arg{[2..Length(arg)]};
        Constraint(Length(L) >= 1); Constraint(ForAll(L, IsSPL));
#        if Length(L) = 1 then return L[1]; fi;
        dims := L[1].dims();
        if not (IsSymbolic(dims[1]) or IsSymbolic(dims[2])) and
           not ForAll(Drop(L, 1), x -> let(d:=x.dims(),
                   IsSymbolic(d[1]) or IsSymbolic(d[2]) or d = dims))
            then Error("Dimensions of summands do not match"); fi;
        return SPL(WithBases(self, rec( simt_dim := simt_dim, _children := L, dimensions := dims)));
    end,

    from_rChildren := (self, rch) >> ApplyFunc(ObjId(self), [self.simt_dim]::rch).appendAobj(self),

    print := meth(self, i, is)
        local s, first, newline;

        if self._short_print or ForAll(self._children, x->not IsRec(x) or IsSPLSym(x) or IsSPLMat(x)) then
            newline := Ignore;
        else
            newline := self._newline;
        fi;
        first := true;
        Print(self.name, "( ", self.simt_dim, ",");
        for s in self._children do
                if(first) then first:=false;
                else Print(", "); fi;
                newline(i + is);
                When(IsSPL(s) or (IsRec(s) and IsBound(s.print) and NumGenArgs(s.print)=2),
                    s.print(i + is, is),
                    Print(s));
        od;
        newline(i);
        Print(")");
        self.printA();

        if IsBound(self._setDims) then
            Print(".overrideDims(", self._setDims, ")");
        fi;
    end,

    ));

Class(SIMTDirectSum, DirectSum, rec(

    abbrevs := [ arg ->
        [ Filtered(
            Flat(List(Flat(arg),
                s -> When(IsSPL(s) and ObjId(s)=DirectSum, s.children(), s))), x-> IsASIMTDim(x) or x.dimensions<>[0,0]) ]
    ],

    new := meth(self, arg)
        local simt_dim, L;
        Constraint(IsList(arg) and Length(arg) >= 2);
        simt_dim := arg[1];
        L := arg{[2..Length(arg)]};
        if ForAll(L, IsIdentitySPL)
            then return I(Sum(L, Rows));
        else
            return SPL(WithBases( self,
            rec( simt_dim := simt_dim, _children := L,
            dimensions := [ Sum(L, t -> t.dimensions[1]),
                            Sum(L, t -> t.dimensions[2]) ] )));
        fi;
    end,

    from_rChildren := (self, rch) >> ApplyFunc(ObjId(self), [self.simt_dim]::rch).appendAobj(self)

    ));

Class(SIMTVStack, VStack, rec(

    abbrevs := [ (simt_dim, spl) -> When(IsASIMTDim(simt_dim) and ObjId(spl) = VStack, [Flat([simt_dim, spl._children])], [[simt_dim, spl]] )  ]
                :: VStack.abbrevs,

    new := (self, arg) >> let(simt_dim := Checked(IsList(arg), arg[1]), spls := arg{[2..Length(arg)]},
                                Checked(
                                    Length(spls) >= 1,
                                    ForAll(spls, s -> IsSymbolic(Cols(s)) or IsSymbolic(Cols(spls[1])) or Cols(s) <= Cols(spls[1])),
                                    SPL(WithBases(self, rec(simt_dim := simt_dim, _children := spls)).setDims())
                                )
                            ),

    from_rChildren := (self, rch) >> ApplyFunc(ObjId(self), [self.simt_dim]::rch).appendAobj(self)

    ));

Class(SIMTHStack, HStack, rec(

    abbrevs := [ (simt_dim, spl) -> When(IsASIMTDim(simt_dim) and ObjId(spl) = HStack, [Flat([simt_dim, spl._children])], [[simt_dim, spl]] )  ]
                :: HStack.abbrevs,

    new := (self, arg) >> let(simt_dim := Checked(IsList(arg), arg[1]), spls := arg{[2..Length(arg)]},
                                Checked(
                                    IsList(spls),
                                    Length(spls) >= 1,
                                    ForAll(spls, s->Rows(s)<=Rows(spls[1])),
                                    SPL(WithBases(self, rec(simt_dim := simt_dim, _children := spls)).setDims())
                                )
                            ),

    from_rChildren := (self, rch) >> ApplyFunc(ObjId(self), [self.simt_dim]::rch).appendAobj(self)

));

Class(SIMTHStack1, SIMTHStack);

Class(SIMTBaseIterative, BaseIterative, rec(

    abbrevs := [ (simt_dim, spl) -> When(IsASIMTDim(simt_dim) and IsIterative(spl), [simt_dim, spl.var, spl.domain, spl._children[1]], [simt_dim, spl] )  ],

    new := meth(self, simt_dim, var, domain, expr)
        local res;
        Constraint(IsSPL(expr));
        # if domain is an integer (not symbolic) it must be positive
        Constraint(not IsInt(domain) or domain >= 0);
        var.isLoopIndex := true;
        res := SPL(WithBases(self, rec(_children := [expr], simt_dim := simt_dim, var := var, domain := domain)));
        res.dimensions := res.dims();
        return res;
    end,

    print := (self, i, is) >> Print(
        self.name, "(", self.simt_dim, ", ",  self.var, ", ", self.domain, ",\n",
        Blanks(i+is), self._children[1].print(i+is, is), "\n",
        Blanks(i), ")", self.printA(),
        When(IsBound(self._setDims), Print(".overrideDims(", self._setDims, ")"), Print(""))
    )
));

Class(SIMTIterVStack, SIMTBaseIterative, rec(
    dims := self >> [ self.child(1).dimensions[1] * self.domain,
                      self.child(1).dimensions[2] ],
));

Class(SIMTIterHStack, SIMTBaseIterative, rec(
    dims := self >> [ self.child(1).dimensions[1],
                      self.child(1).dimensions[2] * self.domain ],
));

Class(SIMTIterHStack1, SIMTIterHStack);

Class(SIMTIDirSum, SIMTBaseIterative, rec(
    dims := self >> let(d:=self._children[1].dimensions, [d[1]* self.domain, d[2]*self.domain])
));

Class(SIMTGrp, Grp);

ext: gi

max_stars_repo_name: broderickpt/spiral-package-simt

max_issues_repo_name: broderickpt/spiral-package-simt

max_forks_repo_name: broderickpt/spiral-package-simt


text:
# Built-in
a := 2988348162058574136915891421498819466320163312926952423791023078876139;
b := 2351399303373464486466122544523690094744975233415544072992656881240319;
PowerModInt(a, b, 10^40);
1527229998585248450016808958343740453059

# Implementation
PowerModAlt := function(a, n, m)
    local r;
    r := 1;
    while n > 0 do
        if IsOddInt(n) then
            r := RemInt(r*a, m);
        fi;
        n := QuoInt(n, 2);
        a := RemInt(a*a, m);
    od;
    return r;
end;

PowerModAlt(a, b, 10^40);

ext: gap

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


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

Declare(_numBlocks);

#F NumBlocks(<sums>, <blocksize>, <comp>)
#F
#F counts the number of blocks of <blocksize> in <sums>, and <comp> is
#F a boolean that tells the function to only count blocks that actually do
#F some computation. It checks for computation by testing the presence of the
#F Blk() structure.

NumBlocks := (S, blocksize, comp) ->
    _numBlocks(S, blocksize, S.dims()[1], comp);

_numBlocks := function(S, blocksize, size, comp)

    if ObjId(S) = ISum then
        size := size / S.var.range;

        if size <= blocksize then
            return When(comp = false or Collect(S, Blk) <> [], 1, 0);
        fi;
    fi;

    return Sum(S.rChildren(), e -> _numBlocks(e, blocksize, size, comp));
end;

Declare(_numWBlocks);

NumWeightedBlocks := (S, blocksize, comp) ->
    _numWBlocks(S, blocksize, S.dims()[1], S.dims()[1], comp);

_numWBlocks := function(S, blocksize, size, fullsize, comp)
    if ObjId(S) = ISum then
        size := size / S.var.range;

        if size <= blocksize then
            return When(comp = false or Collect(S, Blk) <> [],
                fullsize / size,
                0
            );
        fi;
    fi;

    return Sum(S.rChildren(), e -> _numWBlocks(e, blocksize, size, fullsize, comp));
end;

Declare(_numW2Blocks);

#F NumW2Blocks(<sums>, <blocksize>, <comp>)
#F
#F
#F
NumW2Blocks := (S, blocksize, comp) ->
    _numW2Blocks(S, blocksize, S.dims()[1], S.dims()[1], comp);

_numW2Blocks := function(S, blocksize, size, fullsize, comp)

    if ObjId(S) = ISum then
        size := size / S.var.range;

        # <fullsize / size> gives the number of times a block of size
        # <size> is executed.  <2 * size> is the number of out of
        # block accesses a <size> sized block requires. <(2 * size) *
        # (fullsize / size)> is the total number of out-of-block
        # accesses required by this block.
        #
        # <size * fullsize> is an algebraic simplification.
        #
        # since fullsize is a constant factor, we drop it.
        #
        if size <= blocksize then
            return When(comp = false or Collect(S, Blk) <> [],
                size,
                0
            );
        fi;
    fi;

    return Sum(S.rChildren(), e -> _numW2Blocks(e, blocksize, size, fullsize, comp));
end;


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

# Don't forget to Import(platforms.cellSPU)

_svctSizes := (N, p, v) -> Filtered([1..N], i->ForAll(Factors(i), j->j<=p) and IsInt(i/v^2));
_defaultSizes := (func, default) -> ((arg) -> When(Length(arg) = 1, func(When(IsList(arg[1]), arg[1], [arg[1]])), func(default)));

Class(altivecOpts, rec(
            4x32f := rec(
                dft_sc := rec(
                    small := rec(verify:=true, interleavedComplex := false, stdTTensor := false, globalUnrolling:=10000),
                    medium := rec(tsplRader:=false, tsplBluestein:=false, tsplPFA:=false, oddSizes:=false, interleavedComplex := false),
                ),
                dft_ic := rec(
                    small := rec(verify:=true, interleavedComplex := true, stdTTensor := false, globalUnrolling:=10000),
                    medium := rec(tsplRader:=false, tsplBluestein:=false, tsplPFA:=false, oddSizes:=false, interleavedComplex := true),
                )
            )
));

benchAltiVec := function()
    if LocalConfig.cpuinfo.SIMD().hasAltiVec then
        return rec(
            4x32f := rec(
                1d := rec(
                    dft_sc := rec(
                        small := _defaultSizes(s->doSimdDft(s, altivec_4x32f, altivecOpts.4x32f.dft_sc.small), [ 2..64 ]),
                        medium := _defaultSizes(s->doSimdDft(s, altivec_4x32f, altivecOpts.4x32f.dft_sc.medium), _svctSizes(1024, 16, 4)),
#                        medium := spiral.libgen.doParSimdDft(altivec_4x32f, 1, _svctSizes(1024, 16, 4), false, true, false, false),
#                        large := spiral.libgen.doParSimdDft(altivec_4x32f, 1, List([4..20], i->2^i), false, true, false, false)

                    ),
                    dft_ic := rec(
                        small := _defaultSizes(s->doSimdDft(s, altivec_4x32f, altivecOpts.4x32f.dft_ic.small), [ 2..64 ]),
                        medium := _defaultSizes(s->doSimdDft(s, altivec_4x32f, altivecOpts.4x32f.dft_ic.medium), _svctSizes(1024, 16, 4))
#                        medium := spiral.libgen.doParSimdDft(altivec_4x32f, 1, _svctSizes(1024, 16, 4), false, true, false, true),
#                        large := spiral.libgen.doParSimdDft(altivec_4x32f, 1, List([4..20], i->2^i), false, true, false, true)
                    ),
                    rdft := _defaultSizes(s->doSimdSymDFT(TRDFT, s, altivec_4x32f, rec(verify :=true)), 32*[1..12]),
                    dht  := _defaultSizes(s->doSimdSymDFT(TDHT,  s, altivec_4x32f, rec(verify :=true)), 32*[1..12]),
                    dct2 := _defaultSizes(s->doSimdSymDFT(TDCT2, s, altivec_4x32f, rec(verify :=true)), 32*[1..12]),
                    dct3 := _defaultSizes(s->doSimdSymDFT(TDCT3, s, altivec_4x32f, rec(verify :=true)), 32*[1..12]),
                    dct4 := _defaultSizes(s->doSimdSymDFT(TDCT4, s, altivec_4x32f, rec(verify :=true)), 32*[1..12]),
                    dst2 := _defaultSizes(s->doSimdSymDFT(TDST2, s, altivec_4x32f, rec(verify :=true)), 32*[1..12]),
                    dst3 := _defaultSizes(s->doSimdSymDFT(TDST3, s, altivec_4x32f, rec(verify :=true)), 32*[1..12]),
                    dst4 := _defaultSizes(s->doSimdSymDFT(TDST4, s, altivec_4x32f, rec(verify :=true)), 32*[1..12])
                ),

                2d := rec(
                    dft_ic := rec(
                        medium := _defaultSizes(s->doSimdMddft(s, altivec_4x32f, rec(interleavedComplex := true,
                                    oddSizes := false, svct := true, splitL := false, pushTag := true, flipIxA := false, stdTTensor := true, tsplPFA := false)),
                                    16*List([1..5], i->[i,i])),
                        small := _defaultSizes(s->doSimdMddft(s, altivec_4x32f, rec(verify:=true, interleavedComplex := true, globalUnrolling:=10000,
                                    tsplPFA := false, pushTag:= false, oddSizes := true, svct := true, splitL := false)),
                                    List([2..16], i->[i,i]))
                    ),
                    dft_sc := rec(
                        medium := _defaultSizes(s->doSimdMddft(s, altivec_4x32f, rec(interleavedComplex := false,
                                    oddSizes := false, svct := true, splitL := false, pushTag := true, flipIxA := false, stdTTensor := true, tsplPFA := false)),
                                    16*List([1..5], i->[i,i])),
                        small := _defaultSizes(s->doSimdMddft(s, altivec_4x32f, rec(verify:=true, interleavedComplex := false, globalUnrolling:=10000,
                                    tsplPFA := false, pushTag:= false, oddSizes := true, svct := true, splitL := false)),
                                    List([2..16], i->[i,i]))
                    ),
                    dct2 := _defaultSizes(s->doSimdSymMDDFT(DCT2, s, altivec_4x32f, rec(verify := true)), [4, 8, 12, 16]),
                    dct3 := _defaultSizes(s->doSimdSymMDDFT(DCT3, s, altivec_4x32f, rec(verify := true)), [4, 8, 12, 16]),
                    dct4 := _defaultSizes(s->doSimdSymMDDFT(DCT4, s, altivec_4x32f, rec(verify := true)), [4, 8, 12, 16]),
                    dst2 := _defaultSizes(s->doSimdSymMDDFT(DST2, s, altivec_4x32f, rec(verify := true)), [4, 8, 12, 16]),
                    dst3 := _defaultSizes(s->doSimdSymMDDFT(DST3, s, altivec_4x32f, rec(verify := true)), [4, 8, 12, 16]),
                    dst4 := _defaultSizes(s->doSimdSymMDDFT(DST4, s, altivec_4x32f, rec(verify := true)), [4, 8, 12, 16])
                )
            ),
        );
    else
        return false;
    fi;
end;


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

Import(simt);

Class(FFTXCUDAOpts, FFTXOpts, simt.TitanVDefaults, rec(
    tags := [],
    operations := rec(Print := s -> Print("<FFTX CUDA options record>")),
    max_threads := 1024
));

cudaOpts := function(arg)
    local opts;
    opts := Copy(FFTXCUDAOpts);
    opts.breakdownRules.Circulant := [Circulant_PRDFT_FDataNT];
    opts.breakdownRules.PRDFT := List([PRDFT1_Base1, PRDFT1_Base2, PRDFT1_CT, PRDFT1_PF, PRDFT_PD, PRDFT_Rader], _noT);
    opts.breakdownRules.IPRDFT := List([ IPRDFT1_Base1, IPRDFT1_Base2, IPRDFT1_CT, IPRDFT_PD, IPRDFT_Rader ], _noT);
    opts.breakdownRules.PRDFT3 := List([ PRDFT3_Base1, PRDFT3_Base2, PRDFT3_CT ], _noT);
    return opts;

    return opts;
end;

Declare(ParseOptsCUDA);

Class(FFTXCUDADefaultConf, rec(
    getOpts := (self, t) >> ParseOptsCUDA(self, t),
    operations := rec(Print := s -> Print("<FFTX CUDA Default Configuration>")),
    useCUDA := true
));

Class(FFTXCUDADeviceDefaultConf, rec(
    getOpts := (self, t) >> ParseOptsCUDA(self, t),
    operations := rec(Print := s -> Print("<FFTX CUDA Device Default Configuration>")),
    useCUDADevice := true
));

cudaConf := rec(
    defaultName := "defaultCUDAConf",
    defaultOpts := (arg) >> FFTXCUDADefaultConf,
    devFunc := true,
    confHandler := cudaOpts
);

fftx.FFTXGlobals.registerConf(cudaConf);

getTargetOS := function()
    local tgt;

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

#--
Class(FFTXCUDADeviceOpts, FFTXOpts, simt.TitanVDefaults, rec(
    tags := [],
    devFunc := true,
    target := rec ( name := getTargetOS() ),
    operations := rec(Print := s -> Print("<FFTX CUDA Device options record>"))
));

cudaDeviceOpts := function(arg) # specific to WarpX size 100...
    local opts;
    opts := Copy(FFTXCUDADeviceOpts);
    opts.breakdownRules.Circulant := [Circulant_PRDFT_FDataNT];
    opts.breakdownRules.PRDFT := List([PRDFT1_Base1, PRDFT1_Base2, CopyFields(PRDFT1_CT,
            rec(allChildren := P ->Filtered(PRDFT1_CT.allChildren(P), i->When(P[1] = 100, Cols(i[1]) = 4, true)))),
        PRDFT_PD], _noT);
    opts.breakdownRules.IPRDFT := List([ IPRDFT1_Base1, IPRDFT1_Base2, IPRDFT1_CT, IPRDFT_PD ], _noT);
    opts.breakdownRules.PRDFT3 := List([ PRDFT3_Base1, PRDFT3_Base2, PRDFT3_CT ], _noT);
    opts.breakdownRules.DFT := [ DFT_Base,
        CopyFields(DFT_CT, rec(children := nt ->Filtered(DFT_CT.children(nt), i->When(nt.params[1] = 100, Cols(i[1]) = 4, true)))),
        DFT_PD ];
    opts.breakdownRules.TTensorInd := [dsA_base, L_dsA_L_base, dsA_L_base, L_dsA_base];
    return opts;
end;

cudaDeviceConf := rec(
    defaultName := "defaultCUDADeviceConf",
    defaultOpts := (arg) >> FFTXCUDADeviceDefaultConf,
    confHandler := cudaDeviceOpts
);

fftx.FFTXGlobals.registerConf(cudaDeviceConf);

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

    # all dimensions need to be inthis array for the high perf MDDFT conf to kick in for now
    # size 320 is problematic at this point and needs attention. Need support for 3 stages to work first

    _HPCSupportedSizesCUDA := [80, 96, 100, 224, 320];
    _thold := 16;

    if IsBound(conf.useCUDADevice) then
        # detect real MD convolution
        _tt := Collect(t, RCDiag)::Collect(t, MDPRDFT)::Collect(t, IMDPRDFT)::Collect(t, TTensorI);
        if Length(_tt) = 4 then
            _conf := FFTXGlobals.confWarpXCUDADevice();
            _opts := FFTXGlobals.getOpts(_conf);
            return _opts;
        fi;

        # detect batch of DFT/PRDFT/MDDFT/MDPRDFT
        if ((Length(Collect(t, TTensorInd)) >= 1) or (Length(Collect(t, TTensorI)) >= 1)) and
            ((Length(Collect(t, DFT)) = 1) or (Length(Collect(t, PRDFT)) = 1) or (Length(Collect(t, IPRDFT)) = 1) or
              (Length(Collect(t, MDDFT)) >= 1) or (Length(Collect(t, MDPRDFT)) >= 1) or (Length(Collect(t, IMDPRDFT)) >= 1)) then
            _conf := FFTXGlobals.confBatchFFTCUDADevice();
            _opts := FFTXGlobals.getOpts(_conf);
            return _opts;
        fi;

        # detect 3D DFT
        _tt := Collect(t, MDDFT)::Collect(t, MDPRDFT)::Collect(t, IMDPRDFT);
        if Length(_tt) = 1 and Length(_tt[1].params[1]) = 3 then
            _conf := FFTXGlobals.confFFTCUDADevice();
            _opts := FFTXGlobals.getOpts(_conf);

            # opts for high performance CUDA cuFFT
            if Length(Filtered(_tt, i -> ObjId(i) = MDDFT)) > 0 and ForAll(_tt[1].params[1], i-> i in _HPCSupportedSizesCUDA) then
                _opts.breakdownRules.MDDFT := [fftx.platforms.cuda.MDDFT_tSPL_Pease_SIMT];
                _opts.breakdownRules.TTwiddle := [ TTwiddle_Tw1 ];
                _opts.tags := [ASIMTKernelFlag(ASIMTGridDimX), ASIMTBlockDimY, ASIMTBlockDimX];

                _opts.globalUnrolling := 2*_thold + 1;

                _opts.breakdownRules.TTensorI := [CopyFields(IxA_L_split, rec(switch := true)), fftx.platforms.cuda.L_IxA_SIMT]::_opts.breakdownRules.TTensorI;
                _opts.breakdownRules.DFT := [CopyFields(DFT_tSPL_CT, rec(switch := true,
                    filter := e-> When(e[1]*e[2] <= _thold^2, e[1] <= _thold and e[2] <= _thold, e[1] <= _thold and e[2] >= _thold)))]::_opts.breakdownRules.DFT;

                _opts.unparser.simt_synccluster := _opts.unparser.simt_syncblock;
                _opts.postProcessSums := (s, opts) -> let(s1 := ApplyStrategy(s, [ MergedRuleSet(RulesFuncSimp, RulesSums, RulesSIMTFission) ], BUA, opts),
                    FixUpCUDASigmaSPL_3Stage(s1, opts));

                _opts.operations.Print := s -> Print("<FFTX CUDA HPC MDDFT options record>");

            fi;

            return _opts;
        fi;

        # promote with default conf rules
        tt := _promote1(Copy(t));

        if ObjId(tt) = TFCall then
            _tt := tt.params[1];
            # check for convolution
            if (ObjId(_tt) in [MDRConv, MDRConvR]) or ((ObjId(_tt) = TTensorI) and (ObjId(_tt.params[1]) in [MDRConv, MDRConvR])) then
                _conf := FFTXGlobals.confMDRConvCUDADevice();
                _opts := FFTXGlobals.getOpts(_conf);
                return _opts;
            fi;
            # check for Hockney. This is for N=130
            if ObjId(_tt) = IOPrunedMDRConv  and _tt.params[1] = [130,130,130] then
                _conf := FFTXGlobals.confHockneyMlcCUDADevice();
                _opts := FFTXGlobals.getOpts(_conf);
                return _opts;
            fi;
            # check for general Hockney.
            if ObjId(_tt) = IOPrunedMDRConv then
                _conf := FFTXGlobals.confMDRConvCUDADevice();
                _opts := FFTXGlobals.getOpts(_conf);
                _opts.tags := [ASIMTKernelFlag(ASIMTGridDimY), ASIMTGridDimX, ASIMTBlockDimZ];
                return _opts;
            fi;
        fi;

        # check for WarpX
        _conf := FFTXGlobals.confWarpXCUDADevice();
        _opts := FFTXGlobals.getOpts(_conf);
        tt := _opts.preProcess(Copy(t));
        if ObjId(tt) = TFCall and ObjId(tt.params[1]) = TCompose then
            _tt := tt.params[1].params[1];
            # detect promoted WarpX
            if IsList(_tt) and Length(_tt) = 3 and List(_tt, ObjId) = [ TNoDiagPullinRight, TRC, TNoDiagPullinLeft ] then
                return _opts;
            fi;
        fi;
        # we are doing nothing special
        return FFTXGlobals.getOpts(conf);
    fi;
    if IsBound(conf.useCUDA) then
        return FFTXGlobals.getOpts(conf);
    fi;

    # Here we have to handle GPU configs
    Error("Don't know how to derive opts!\n");
end;

max_stars_repo_name: broderickpt/spiral-package-fftx

max_issues_repo_name: broderickpt/spiral-package-fftx

max_forks_repo_name: broderickpt/spiral-package-fftx


text:
# Determine images of the homomorphism between to a permutation group
# and a matrix groups both isomorphic to S_4.
G := Group((1, 2, 3, 4), (2, 3, 4));
H := Group([[0, -1, 0], [1, 0, 0], [0, 0, 1]], [[0, 1, 0], [0, 0, 1], [1, 0, 0]]);

hom := GroupHomomorphismByImages(G, H);

for g in G do
    l := [1..4];
    Apply(l, i -> String((i^g) - 1));
    w := Concatenation(l);
    image := Image(hom, g);
    Print(w, " ", image, "\n");
od;

max_stars_repo_name: fifth-postulate/packing-puzzle

max_issues_repo_name: fifth-postulate/packing-puzzle

max_forks_repo_name: fifth-postulate/packing-puzzle


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

NewRulesFor(DFT, rec(
    #F DFT_Base: DFT_2 = F_2
    #F
    DFT_tSPL_VBase := rec(
        switch := false,
        applicable := (t) -> Rows(t) = 2 and t.isTag(1,AVecReg),
        children := t -> [[ TTensorI(DFT(2), 1, AVec, AVec).withTags(t.getTags()) ]],
        apply := (t, C, Nonterms) -> C[1]
    ),
    #F DFT_Base: DFT_2 = F_2 for complex v/2-way vectorization
    #F
    DFT_tSPL_CxVBase := rec(
        switch := false,
        applicable := (t) -> Rows(t) = 2 and t.isTag(1,AVecRegCx),
        apply := (t, C, Nonterms) -> let(
	    v := t.getTags()[1].v,
	    Cond(v = 1, VTensor(F(2), 1),
		 # else
		 let(zeros    := Replicate(v-2, 0.0),
		     one_one  := TVect(TReal, v).value( [1.0,  1.0] :: zeros ),
		     one_mone := TVect(TReal, v).value( [1.0, -1.0] :: zeros ),
		     blk      := VBlk([[ one_mone, one_one ]], v),
		     scat     := _VVStack([VTensor(I(1), v), VIxJ2(v)], v),
		     IxVScat_pc(1,2,2,0,v) * blk * scat * IxVGath_pc(1,2,2,0,v)))
        )
    ),

    #F DFT_Base: DFT_4 for complex v/2-way vectorization
    #F Tensor(F(2),I(4)) * RC(Diag(Tw1(4,2,-1))) * Tensor(I(2), F(2), I(2)) * Tensor(L(4,2),I(2))
    DFT_tSPL_CxVBase2 := rec(
        switch := false,
        applicable := (t) -> t.params[1] = 4 and t.firstTagIs(AVecRegCx) and t.firstTag().v>=4,
        children := nt -> [[TL(4, 2, 1, 1).withTags(nt.getTags()).setWrap(VWrapId)]],

        apply := (t, C, Nonterms) -> let(
	    v  := t.firstTag().v,
            vt := TVect(TReal, v),

            IxVScat_pc(1,4,4,0,v) *
            VBlk( [[ vt.value([1, 1, -1,-1]), vt.one() ]], v) *
            _VVStack([VTensor(I(1), v), VJxI(2, v)], v) *
            VDiag(Tw1(4,2,t.params[2]), v) *
            VBlk( [[ vt.value([1, -1, 1,-1]), vt.one() ]], v) *
            _VVStack([VTensor(I(1), v), VIxJ2(v)], v) *
            C[1] *
            IxVGath_pc(1,4,4,0,v)
        )
    ),
    #F DFT_tSPL_CT_oddvloop
    #F DFT(mn, k)->Tensor(t2 * RI(m, mu), RI(n,nu)) * (Tensor(RI(mu,m), RI(nu,n)) * Diag(twf) * Tensor(RI(m,mu), RI(n,nu))) * L(mu*nu, mu) * Tensor(RI(nu,n) * t1, RI(mu, m))
    DFT_tSPL_CT_oddvloop := rec(
        switch  := false,
        maxSize := false,
        filter  := e->true,
        applicable := (self, nt) >> nt.params[1] > 2
            and (self.maxSize = false or nt.params[1] <= self.maxSize)
            and not IsPrime(nt.params[1])
            and (nt.firstTagIs(AVecReg) or nt.firstTagIs(AVecRegCx)),

        children := (self, nt) >> let(
            pv := nt.getTags(),
            v := pv[1].v,
            isa := pv[1].isa,
            _wrp := _nt -> When(Length(pv)=1, _nt.setWrap(VWrap(isa)), _nt.withTags(Drop(pv, 1)).setWrap(VWrap(isa))),
            Map2(Filtered(Filtered(DivisorPairs(nt.params[1]), a->IsOddInt(a[1])), self.filter),
                (m,n) -> [ _wrp(DFT(m, nt.params[2] mod m)),
                            _wrp(DFT(n, nt.params[2] mod n)),
                            TL(_roundup(n,v).v*v,v).withTags(pv).setWrap(VWrapId)
                         ])),

        apply := (t, C, Nonterms) -> let(
  	        v := t.firstTag().v,
            lprm := C[3],
            m := Rows(C[1]),
            n := Rows(C[2]),
            N := m * n,
            k := t.params[2],
            dftm := C[1],
            dftn := C[2],
            mu := _roundup(m,v).v,
            md := _rounddown(m,v).v,
            nu := _roundup(n,v).v,
            nd := _rounddown(n,v).v,
            f1 := fCompose(Tw1(N, n, k), fTensor(fId(m), fAdd(n,nd,0))),
            f2 := When(mu = m,
                fCompose(Tw1(N, n, k), fTensor(fId(m), fAdd(n,n-nd,nd))),
                fStretch(fCompose(Tw1(N, n, k), fTensor(fId(m), fAdd(n,n-nd,nd))), v, mu-m)),
            i1 := Ind(md/v),
            i2 := Ind(nd/v),
            #---------------------
            stage1 := ApplyFunc(SUM, Concat([
                ISum(i1,
                    VScat(fTensor(fAdd(mu/v,1,i1), fId(nu)),v) *
                    lprm * VScat_zero(nu, n, v) * VTensor(dftn,v) *
                    IxVGath_pc(n, m, v, v*i1, v)
                )],
                Cond(IsOddInt(m),
                    [VScat(fTensor(fAdd(mu/v,1,md/v), fId(nu)),v) *
                    lprm * VScat_zero(nu, n, v) * VTensor(dftn,v) *
                    IxVGath_pc(n, m, m-md, md, v)],
                    [])
                )
            ),
            #---------------------
            stage2 := ApplyFunc(SUM, Concat([
                ISum(i2,
                    IxVScat_pc(m, n, v, v*i2, v) *
                    VTensor(dftm, v) * VDiag(fPrecompute(fCompose(f1, fTensor(fId(m), fBase(i2), fId(v)))), v) * VGath(fAdd(mu, m, 0), v) *
                    VGath(fTensor(fId(mu), fAdd(nu/v, 1,i2)),v)
                )],
                Cond(IsOddInt(n),
                    [IxVScat_pc(m, n, n-nd, nd, v) *
                    VTensor(dftm, v) * VDiag(f2, v) * VGath(fAdd(mu, m, 0), v) *
                    VGath(fTensor(fId(mu), fAdd(nu/v, 1,nd/v)),v)],
                    [])
                )
            ),
#            #---------------------
#            stage1 := SUM(
#                ISum(i1,
#                    VScat(fTensor(fAdd(mu/v,1,i1), fId(nu)),v) *
#                    lprm * VScat_zero(nu, n, v) * VTensor(dftn,v) *
#                    IxVGath_pc(n, m, v, v*i1, v)
#                ),
#                VScat(fTensor(fAdd(mu/v,1,md/v), fId(nu)),v) *
#                lprm * VScat_zero(nu, n, v) * VTensor(dftn,v) *
#                IxVGath_pc(n, m, m-md, md, v)
#            ),
#            #---------------------
#            stage2 := SUM(
#                ISum(i2,
#                    IxVScat_pc(m, n, v, v*i2, v) *
#                    VTensor(dftm, v) * VDiag(fPrecompute(fCompose(f1, fTensor(fId(m), fBase(i2), fId(v)))), v) * VGath(fAdd(mu, m, 0), v) *
#                    VGath(fTensor(fId(mu), fAdd(nu/v, 1,i2)),v)
#                ),
#                IxVScat_pc(m, n, n-nd, nd, v) *
#                VTensor(dftm, v) * VDiag(f2, v) * VGath(fAdd(mu, m, 0), v) *
#                VGath(fTensor(fId(mu), fAdd(nu/v, 1,nd/v)),v)
#            ),
            #---------------------
            stage2 * stage1
        )
    ),

#    #F DFT_tSPL_CT_evenoddvloop
#    #F DFT(mn, k)->Tensor(t2 * RI(m, mu), RI(n,nu)) * (Tensor(RI(mu,m), RI(nu,n)) * Diag(twf) * Tensor(RI(m,mu), RI(n,nu))) * L(mu*nu, mu) * Tensor(RI(nu,n) * t1, RI(mu, m))
#    DFT_tSPL_CT_evenoddvloop := rec(
#        switch  := false,
#        maxSize := false,
#        filter  := e->true,
#        applicable := (self, nt) >> nt.params[1] > 2
#            and (self.maxSize = false or nt.params[1] <= self.maxSize)
#            and not IsPrime(nt.params[1])
#            and IsOddInt(nt.params[1]/2)
#            and IsEvenInt(nt.params[1])
#            and (nt.firstTagIs(AVecReg) or nt.firstTagIs(AVecRegCx)),
#
#        children := (self, nt) >> let(
#            pv := nt.getTags(),
#            v := pv[1].v,
#            isa := pv[1].isa,
#            _wrp := _nt -> When(Length(pv)=1, _nt.setWrap(VWrap(isa)), _nt.withTags(Drop(pv, 1)).setWrap(VWrap(isa))),
#            Map2(Filtered(DivisorPairs(nt.params[1]), self.filter),
#                (m,n) -> [ _wrp(DFT(m, nt.params[2] mod m)),
#                            _wrp(DFT(n, nt.params[2] mod n)),
#                            TL(_roundup(n,v).v*v,v).withTags(pv).setWrap(VWrapId)
#                         ])),
#
#        apply := (t, C, Nonterms) -> let(
#  	        v := t.firstTag().v,
#            lprm := C[3],
#            m := Rows(C[1]),
#            n := Rows(C[2]),
#            N := m * n,
#            k := t.params[2],
#            dftm := C[1],
#            dftn := C[2],
#            mu := _roundup(m,2).v,
#            md := _rounddown(m,2).v,
#            nu := _roundup(n,2).v,
#            nd := _rounddown(n,2).v,
#            f1 := fCompose(Tw1(N, n, k), fTensor(fId(m), fAdd(n,nd,0))),
#            f2 := fStretch(fCompose(Tw1(N, n, k), fTensor(fId(m), fAdd(n,n-nd,nd))), v, mu-m),
#            i1 := Ind(md/v),
#            i2 := Ind(nd/v),
#            #---------------------
#            stage1 := ApplyFunc(SUM, Concat([
#                ISum(i1,
#                    VScat(fTensor(fAdd(mu/v,1,i1), fId(nu)),v) *
#                    lprm * VScat_zero(nu, n, v) * VTensor(dftn,v) *
#                    IxVGath_pc(n, m, v, v*i1, v)
#                )],
#                Cond(IsOddInt(m),
#                    [VScat(fTensor(fAdd(mu/v,1,md/v), fId(nu)),v) *
#                    lprm * VScat_zero(nu, n, v) * VTensor(dftn,v) *
#                    IxVGath_pc(n, m, m-md, md, v)],
#                    [])
#                )
#            ),
#            #---------------------
#            stage2 := ApplyFunc(SUM, Concat([
#                ISum(i2,
#                    IxVScat_pc(m, n, v, v*i2, v) *
#                    VTensor(dftm, v) * VDiag(fPrecompute(fCompose(f1, fTensor(fId(m), fBase(i2), fId(v)))), v) * VGath(fAdd(mu, m, 0), v) *
#                    VGath(fTensor(fId(mu), fAdd(nu/v, 1,i2)),v)
#                )],
#                Cond(IsOddInt(n),
#                    [IxVScat_pc(m, n, n-nd, nd, v) *
#                    VTensor(dftm, v) * VDiag(f2, v) * VGath(fAdd(mu, m, 0), v) *
#                    VGath(fTensor(fId(mu), fAdd(nu/v, 1,nd/v)),v)],
#                    [])
#                )
#            ),
#            #---------------------
#            stage2 * stage1
#        )
#    )
));


text:
#############################################################################
####
##
#W  anupqxdesc.gi              ANUPQ package                    Werner Nickel
#W                                                                Greg Gamble
##
##  Installs functions to do recursive development of a descendants tree.
##  If ANUPQ is loaded from XGAP the development is seen graphically.
##
#Y  Copyright (C) 2001  Lehrstuhl D fuer Mathematik,  RWTH Aachen,  Germany
##

#############################################################################
##
#F  PqDescendantsTreeCoclassOne([<i>]) . . . generate a coclass one des. tree
##
##  for the <i>th  or  default  interactive  {\ANUPQ}  process,  generates  a
##  descendant tree for the  group  of  the  process  (which  must  be  a  pc
##  $p$-group) consisting of descendants of $p$-coclass 1  and  extending  to
##  the class determined by the option `TreeDepth' (or 6  if  the  option  is
##  omitted). In an  {\XGAP}  session,  a  graphical  representation  of  the
##  descendants tree appears  in  a  separate  window.  Subsequent  calls  to
##  `PqDescendantsTreeCoclassOne' for the same process may be used to  extend
##  the descendant tree from the last descendant  computed  that  itself  has
##  more than one descendant. `PqDescendantsTreeCoclassOne' also accepts  the
##  options  `CapableDescendants'  (or  `AllDescendants')  and  any   options
##  accepted     by     the     interactive     `PqDescendants'      function
##  (see~"PqDescendants!interactive").
##
##  *Notes*
##
##  `PqDescendantsTreeCoclassOne'    first    calls    `PqDescendants'.    If
##  `PqDescendants' has already been called for  the  process,  the  previous
##  value computed is used and a warning is `Info'-ed at `InfoANUPQ" level 1.
##
##  As each descendant is processed its unique  label  defined  by  the  `pq'
##  program and number of descendants is `Info'-ed at `InfoANUPQ' level 1.
##
##  `PqDescendantsTreeCoclassOne' is an  ``experimental''  function  that  is
##  included to demonstrate the sort of things that  are  possible  with  the
##  $p$-group generation machinery.
##
InstallGlobalFunction( PqDescendantsTreeCoclassOne, function( arg )
    local   datarec,  title,  des,  node;

    datarec := ANUPQ_ARG_CHK("PqDescendantsTreeCoclassOne", arg);
    if datarec.procId = 0 then
        Error("non-interactive `PqDescendantsTreeCoclassOne' is not ",
              "currently supported\n");
    fi;
    if IsBound(datarec.treepos) then
        PQX_RECURSE_DESCENDANTS( datarec,
                                 datarec.treepos.class,
                                 datarec.treepos.node,
                                 datarec.treepos.ndes );
        return;
    fi;

    des := PqDescendants( datarec.procId : StepSize := 1 );
    if IsPackageMarkedForLoading("xgap","") then
        title := Concatenation( "Descendants Tree, p=",
                     String( PrimePGroup(datarec.group) ),
                     ", order: ",
                     String( Size(datarec.group) ),
                     ", class <= ",
                     String( VALUE_PQ_OPTION("TreeDepth", 6, datarec) ) );
        if VALUE_PQ_OPTION("CapableDescendants",
                           not VALUE_PQ_OPTION("AllDescendants", true, datarec),
                           datarec) then
            Append( title, ", capable descendants" );
        else
            Append( title, ", all descendants" );
        fi;
        datarec.xgapsheet := GraphicSheet( title, 800, 700 );
        datarec.nextX     := [0, 0];
        node := PQX_PLACE_NEXT_NODE( datarec, 1 );
    else
        node := 0;
    fi;

    PQX_RECURSE_DESCENDANTS( datarec, 2, node, Length(des) );

    if datarec.calltype = "non-interactive" then
        PQ_COMPLETE_NONINTERACTIVE_FUNC_CALL(datarec);
        if IsBound( datarec.setupfile ) then
            return true;
        fi;
    fi;
end);

#############################################################################
##
#F  PQX_PLACE_NEXT_NODE( <datarec>, <class> ) . place a node on an XGAP sheet
##
##  places a node for the current descendant of class <class> on the  {\XGAP}
##  sheet `<datarec>.xgapsheet'.
##
InstallGlobalFunction( PQX_PLACE_NEXT_NODE, function( datarec, class )
    local   y,  x;

    y := 40 * class;
    datarec.nextX[class] := datarec.nextX[class] + 16;
    x := datarec.nextX[class];

    Disc( datarec.xgapsheet, x, y, 6 );
    return [x,y];
end);

#############################################################################
##
#F  PQX_MAKE_CONNECTION( <datarec>, <a>, <b> ) . .  join two XGAP sheet nodes
##
##  joins the nodes <a> and <b> on the  {\XGAP}  sheet  `<datarec>.xgapsheet'
##  with a straight line, where <a> and <b> are each `[<x>, <y>]'  coordinate
##  pairs (lists) of integers.
##
InstallGlobalFunction( PQX_MAKE_CONNECTION, function( datarec, a, b )

    Line( datarec.xgapsheet, a[1], a[2], b[1]-a[1], b[2]-a[2] );
end);

#############################################################################
##
#F  PQX_RECURSE_DESCENDANTS(<datarec>,<class>,<parent>,<n>)  extend des. tree
##
##  extends a descendant tree of coclass 1 descendants from the current  `pq'
##  descendant  that  has  <n>  descendants  from  class  <class>  to   class
##  determined by the option `TreeDepth' (or 6, by  default)  from  the  node
##  <parent> which is an `[<x>, <y>]' coordinate pair (list) of integers.
##
InstallGlobalFunction( PQX_RECURSE_DESCENDANTS,
function( datarec, class, parent, n)
local   i,  node,  nr;

    if class > VALUE_PQ_OPTION("TreeDepth", 6, datarec) then
        datarec.treepos := rec(class:=class, node:=parent, ndes:=n);
        return;
    fi;

    if parent <> 0 then
        datarec.nextX[class] := datarec.nextX[class] + 14;
    fi;
    for i in [1..n] do
        PQ_PG_RESTORE_GROUP( datarec, class, i );
        PQ_PG_EXTEND_AUTOMORPHISMS( datarec );
        nr := PqPGConstructDescendants( datarec.procId : StepSize := 1 );

        if parent <> 0 and (not datarec.CapableDescendants or nr > 0) then
            # Place a node on the Graphic Sheet and connect it with its
            # parent node.
            node := PQX_PLACE_NEXT_NODE( datarec, class );
            PQX_MAKE_CONNECTION( datarec, node, parent );
        else
            node := 0;
        fi;

        Info(InfoANUPQ, 1, "Number of descendants of group ", datarec.gpnum,
                           " to class ", class, ": ", nr);
        if nr > 0 then
            if IsBound(datarec.nextX) and
               not IsBound(datarec.nextX[class + 1]) then
                datarec.nextX[class + 1] := 0;
            fi;
            PQX_RECURSE_DESCENDANTS( datarec, class+1, node, nr );
        fi;
    od;
    return;
end);

#E  anupqxdesc.gi . . . . . . . . . . . . . . . . . . . . . . . . . ends here

max_stars_repo_name: gap-system/anupq

max_issues_repo_name: gap-system/anupq

max_forks_repo_name: gap-system/anupq


text:
ReadByLines := function(name)
	local file, line, count;
	file := InputTextFile(name);
	count := 0;
	while true do
		line := ReadLine(file);
		if line = fail then
			break;
		fi;
		count := count + 1;
	od;
	CloseStream(file);
	return count;
end;

# With [http://www.ibiblio.org/pub/docs/misc/amnesty.txt amnesty.txt]
ReadByLines("amnesty.txt");
# 384


text:
leftOf  :=function(setA, vA, setB, vB)
local i;
for i in [1..4] do
   if ( setA[i] = vA) and  (setB[i+1] = vB) then return true ;fi;
od;
 return false;
end;

nextTo  :=function(setA, vA, setB, vB)
local i;
for i in [1..4] do
   if ( setA[i] = vA) and  (setB[i+1] = vB) then return true ;fi;
   if ( setB[i] = vB) and  (setA[i+1] = vA) then return true ;fi;
od;
return false;
end;

requires := function(setA, vA, setB, vB)
    local i;
      for i in [1..5] do
        if ( setA[i] = vA) and  (setB[i] = vB) then return true ;fi;
      od;
 return false;
end;

pcolors :=PermutationsList(["white" ,"yellow" ,"blue" ,"red" ,"green"]);
pcigars :=PermutationsList(["blends", "pall_mall", "prince", "bluemasters", "dunhill"]);
pnats:=PermutationsList(["german", "swedish", "british", "norwegian", "danish"]);
pdrinks :=PermutationsList(["beer", "water", "tea", "milk", "coffee"]);
ppets  :=PermutationsList(["birds", "cats", "horses", "fish", "dogs"]);

for colors in pcolors do
if not (leftOf(colors,"green",colors,"white")) then continue ;fi;
for nats in pnats do
if not (requires(nats,"british",colors,"red")) then  continue ;fi;
if not (nats[1]="norwegian") then continue ;fi;
if not (nextTo(nats,"norwegian",colors,"blue")) then continue ;fi;
for pets in ppets do
if not (requires(nats,"swedish",pets,"dogs")) then  continue ;fi;
for drinks in pdrinks do
if not (drinks[3]="milk") then continue ;fi;
if not (requires(colors,"green",drinks,"coffee")) then continue ;fi;
if not (requires(nats,"danish",drinks,"tea")) then  continue ;fi;
for cigars in pcigars do
if not (nextTo(pets,"horses",cigars,"dunhill")) then continue ;fi;
if not (requires(cigars,"pall_mall",pets,"birds")) then  continue ;fi;
if not (nextTo(cigars,"blends",drinks,"water")) then  continue ;fi;
if not (nextTo(cigars,"blends",pets,"cats")) then  continue ;fi;
if not (requires(nats,"german",cigars,"prince")) then  continue ;fi;
if not (requires(colors,"yellow",cigars,"dunhill")) then continue ;fi;
if not (requires(cigars,"bluemasters",drinks,"beer")) then  continue ;fi;
Print(colors,"\n");
Print(nats,"\n");
Print(drinks,"\n");
Print(pets,"\n");
Print(cigars,"\n");
od;od;od;od;od;


text:
Class(SP, SumsBase, BaseMat, rec(
    _short_print := true,
	rChildren := self >> [self.sr, self.list1, self.list2],
	rSetChild := rSetChildFields("sr","list1", "list2"),
	dims := self >> [self.list1.dimensions[2], self.list2.dimensions[2]],

	new := (self, sr, list1, list2) >> SPL(WithBases(self,
		rec(dimensions := [list1.dimensions[2], list2.dimensions[2]],
		sr := Checked(IsSemiring(sr), sr),
		list1 := list1,
		list2 := list2))),

));

Class(SP2, SumsBase, BaseMat, rec(
    _short_print := true,
	rChildren := self >> [self.list1, self.list2],
	rSetChild := rSetChildFields("list1", "list2"),

	new := (self, list1, list2) >> SPL(WithBases(self,
		rec(
		list1 := list1,
		list2 := list2))),

	intersect := meth(list1, list2)
		local i, x, result;
		x := 1;
		i := 1;
		result := 0;
		while i < Length(list1) and x < Length(list2) do
			if list1[i] < list2[x] then
				i := i + 1;
			elif list2[x] < list1[i] then
				x := x + 1;
			else
				result := result + 1;
				i := i +1;
				x := x +1;
			fi;
		od;
		return result;
	end,
));

Class(Intersect, BaseOperation, rec(
  __call__ := (self, sr, x_n, s1, s2, range) >> SPL(WithBases(self, rec(
    sr := Checked(IsSemiring(sr), sr),
    start1 := s1,
    start2 := s2,
    range := range,
    x_n := x_n,
    _children := [range]
  ))),
  dims := self >> [self.x_n, V(1)],
  from_rChildren := (self, rch) >> CopyFields(self, rec(_children := rch)),
  doNotMarkBB := true,
  print := (self, i, is) >> Print(
        self.name, "(", self.sr, ", ", self.start1, ", ", self.start2, ", ", self.range, ")"
  )
));

Class(ScalarProd, TaggedNonTerminal, rec(
  abbrevs := [ (n) -> [n] ],
  isReal := self >> true,
  dims := self >> [self.params[1], self.params[1]]
));

NewRulesFor(ScalarProd, rec(
  Intersect := rec(
    info := "set intersection across two vectors",
    maxSize := false,
    applicable := nt -> not nt.hasTags(),
    children := (self, nt) >> let(i := Ind(), n := nt.params[1], sr := TSemiring_Arithmetic(n.t),
              [[ Intersect(sr, [i, V(0)], [V(0), V(0)], [i, V(0)], i, i, n)]]),
    apply := (nt, C, ct) -> let(i := C[1]._children[1].rt.node.d1, n := nt.params[1], sr := TSemiring_Arithmetic(n.t),C[1])
  ),
));

DefaultSumsGen.Intersect := (self, o, opts) >> o;

DefaultSumsGen.SP := (self, o, opts) >> o;

DefaultCodegen.SP := meth(self, o, y, x, opts)
	local start, iA, jA, e1, istart, jstart, kstart, kend;
	iA := var.fresh_t("iA", TPtr(x.t));
	jA := var.fresh_t("jA", TPtr(x.t));
	start := var.fresh_t("start", TPtr(x.t));
	e1 := var.fresh_t("e1", TPtr(x.t));
	istart := var.fresh_t("istart", TPtr(x.t));
	jstart := var.fresh_t("jstart", TPtr(x.t));
	kstart := var.fresh_t("kstart", TPtr(x.t));
	kend := var.fresh_t("kend", TPtr(x.t));
	return decl([jstart,start, kstart,kend],chain(
	assign(y, o.sr.zero()),
	assign(jstart, start),
	assign(kstart, add(jA, deref(add(iA, deref(istart)))),
	assign(kend, add(jA, deref(add(iA, add(deref(istart),V(1)))))),
	loopw(logic_and(lt(jstart, e1), lt(kstart,kend)),
	chain(
		if3(lt(deref(jstart), deref(kstart)),
			assign(jstart, add(jstart,V(1))),
			lt(deref(kstart), deref(jstart)),
			assign(kstart, add(kstart, V(1))),
			chain(
				assign(y, o.srt.sum(y, V(1))),
				assign(jstart, add(jstart, V(1))),
				assign(kstart, add(kstart, V(1))))))))));
end;

TCDefaults := CopyFields(SpiralDefaults, rec(
  compileStrategy := GraphIndicesCS,
  X := var("IJ", TPtr(TInt)),
  XType := TPtr(TInt),
  arrayDataModifier := "",
  arrayBufModifier := "",
  Y := var("res", TPtr(TInt)),
  YType := TPtr(TInt),
  isCSR := true,
  includes := ["<include/sparse.h>"],
));


text:
#############################################################################
####
##
#W  anupqios.gi            ANUPQ package                          Greg Gamble
##
##  This file installs core functions used with iostreams.
##
#Y  Copyright (C) 2001  Lehrstuhl D fuer Mathematik,  RWTH Aachen,  Germany
##

#############################################################################
##
#F  PQ_START( <workspace>, <setupfile> ) . . . open a stream for a pq process
##
##  ensures the images file written by the `pq' binary when in  the  Standard
##  Presentation menu is empty, opens an io stream  to  a  `pq'  process  (if
##  <setupfile> is `fail') or a file stream for a setup file (if  <setupfile>
##  is a filename i.e. a string) and returns  a  record  with  fields  `menu'
##  (current menu for the `pq' binary), `opts' (the runtime switches used  by
##  the `pq' process), `workspace' (the value of <workspace> which should  be
##  a positive integer), and `stream' (the io or file stream opened).
##
InstallGlobalFunction(PQ_START, function( workspace, setupfile )
local opts, iorec, topqlogfile;
  PrintTo(ANUPQData.SPimages, ""); #to ensure it's empty
  if setupfile = fail then
    opts := [ "-G" ];
  else
    opts := [ "-i", "-k", "-g" ];
  fi;
  if workspace <> 10000000 then
    Append( opts, [ "-s", String(workspace) ] );
  fi;
  iorec := rec( menu := "SP",
                opts := opts,
                workspace := workspace );
  if setupfile = fail then
    iorec.stream := InputOutputLocalProcess( ANUPQData.tmpdir,
                                             ANUPQData.binary,
                                             opts );
    if iorec.stream = fail then
      Error( "failed to launch child process" );
    fi;
    # menus are flushed at InfoANUPQ level 6, prompts at level 5
    FLUSH_PQ_STREAM_UNTIL(iorec.stream, 6, 5, PQ_READ_NEXT_LINE, IS_PQ_PROMPT);
  else
    iorec.stream := OutputTextFile(setupfile, false);
    iorec.setupfile := setupfile;
    ToPQk(iorec, [], [ "#call pq with flags: '",
                       JoinStringsWithSeparator(opts, " "),
                       "'" ]);
  fi;
  return iorec;
end );

#############################################################################
##
#F  PqStart(<G>,<workspace> : <options>) . Initiate interactive ANUPQ session
#F  PqStart(<G> : <options>)
#F  PqStart(<workspace> : <options>)
#F  PqStart( : <options>)
##
##  activate an iostream for an interactive {\ANUPQ} process (i.e.  `PqStart'
##  starts up a `pq' binary process and opens a {\GAP} iostream  to  ``talk''
##  to that process) and returns an integer <i> that can be used to  identify
##  that process. The argument <G>, if given, should be an *fp group* or  *pc
##  group* that the user  intends  to  manipute  using  interactive  {\ANUPQ}
##  functions. If `PqStart' is given an integer argument <workspace> then the
##  `pq' binary is started up with a workspace (an  integer  array)  of  size
##  <workspace> (i.e. $4 \times <workspace>$ bytes in a 32-bit  environment);
##  otherwise, the `pq' binary sets a default workspace of $10000000$.
##
##  The only <options> currently recognised  by  `PqStart'  are  `Prime'  and
##  `Exponent' (see~"Pq" for details) and if provided  they  are  essentially
##  global for the interactive {\ANUPQ} process, except that any  interactive
##  function interacting with the process and passing new  values  for  these
##  options will over-ride the global values.
##
InstallGlobalFunction(PqStart, function(arg)
local opts, iorec, procId, G, workspace, optname;

  if 2 < Length(arg) then
    Error("at most two arguments expected.\n");
  fi;

  if not IsEmpty(arg) and IsGroup( arg[1] ) then
    G := arg[1];
    if not( IsFpGroup(G) or IsPcGroup(G) ) then
      Error( "argument <G> should be an fp group or a pc group\n" );
    fi;
    arg := arg{[2 .. Length(arg)]};
  fi;

  if not IsEmpty(arg) then
    workspace := arg[1];
    if not IsPosInt(workspace) then
      Error("argument <workspace> should be a positive integer.\n");
    fi;
  else
    workspace := 10000000;
  fi;

  iorec := PQ_START( workspace, fail );
  if IsBound( G ) then
    iorec.group := G;
  fi;
  iorec.calltype := "interactive";
  for optname in ANUPQGlobalOptions do
    VALUE_PQ_OPTION(optname, iorec);
  od;

  procId := Length(ANUPQData.io) + 1;
  iorec.procId := procId;
  ANUPQData.io[ procId ] := iorec;
  return procId;
end);

#############################################################################
##
#F  PqQuit( <i> )  . . . . . . . . . . . . Close an interactive ANUPQ session
#F  PqQuit()
##
##  closes the stream of the <i>th or default  interactive  {\ANUPQ}  process
##  and unbinds its `ANUPQData.io' record.
##
InstallGlobalFunction(PqQuit, function(arg)
local ioIndex;

  ioIndex := ANUPQ_IOINDEX(arg);
  # No need to bother about descending through the menus.
  CloseStream(ANUPQData.io[ioIndex].stream);
  Unbind(ANUPQData.io[ioIndex]);
end);

#############################################################################
##
#F  PqQuitAll() . . . . . . . . . . . .  Close all interactive ANUPQ sessions
##
##  closes the streams of all active interactive {\ANUPQ} process and unbinds
##  their `ANUPQData.io' records.
##
InstallGlobalFunction(PqQuitAll, function()
local ioIndex;

  for ioIndex in [1 .. Length(ANUPQData.io)] do
    if IsBound(ANUPQData.io[ioIndex]) then
      CloseStream(ANUPQData.io[ioIndex].stream);
      Unbind(ANUPQData.io[ioIndex]);
    fi;
  od;
end);

#############################################################################
##
#F  ANUPQ_IOINDEX . . . . the number identifying an interactive ANUPQ session
##
##  returns the index of the record in the `ANUPQData.io' list  corresponding
##  to an interactive {\ANUPQ} session. With  no  argument  the  first  bound
##  index in `ANUPQData.io' is returned. With integer (first)  argument  <i>,
##  <i> is returned if `ANUPQData.io[<i>]' is bound.
##
InstallGlobalFunction(ANUPQ_IOINDEX, function(arglist)
local ioIndex;

  if IsEmpty(arglist) then
    # Find the first bound ioIndex
    ioIndex := 1;
    while not(IsBound(ANUPQData.io[ioIndex])) and
          ioIndex < Length(ANUPQData.io) do
      ioIndex := ioIndex + 1;
    od;
    if IsBound(ANUPQData.io[ioIndex]) then
      return ioIndex;
    else
      Info(InfoANUPQ + InfoWarning, 1,
           "No interactive ANUPQ sessions are currently active");
      return fail;
    fi;
  elif IsBound(ANUPQData.io[ arglist[1] ]) then
    return arglist[1];
  else
    Error("no such interactive ANUPQ session\n");
  fi;
end);

#############################################################################
##
#F  ANUPQ_IOINDEX_ARG_CHK .  Checks ANUPQ_IOINDEX has the right no. of arg'ts
##
InstallGlobalFunction(ANUPQ_IOINDEX_ARG_CHK, function(arglist)
  if Length(arglist) > 1 then
    Info(InfoANUPQ + InfoWarning, 1,
         "Expected 0 or 1 arguments, all but first argument ignored");
  fi;
end);

#############################################################################
##
#F  ANUPQDataRecord([<i>]) . . . . . . . returns the data record of a process
##
InstallGlobalFunction(ANUPQDataRecord, function( arg )
  if not IsEmpty(arg) and arg[1] = 0 and IsBound( ANUPQData.ni ) then
    return ANUPQData.ni;
  else
    return ANUPQData.io[ CallFuncList(PqProcessIndex, arg) ];
  fi;
end);

#############################################################################
##
#F  PqProcessIndex( <i> ) . . . . . . . . . . . User version of ANUPQ_IOINDEX
#F  PqProcessIndex()
##
##  If given (at least) one integer  argument  `PqProcessIndex'  returns  its
##  first argument if it corresponds to  an  active  interactive  process  or
##  raises an error; otherwise, with no arguments,  it  returns  the  default
##  active interactive process. If the user provides more than  one  argument
##  then all arguments other than the  first  argument  are  ignored  (and  a
##  warning is issued to `Info' at `InfoANUPQ' or `InfoWarning' level 1).
##
InstallGlobalFunction(PqProcessIndex, function(arg)
  ANUPQ_IOINDEX_ARG_CHK(arg);
  return ANUPQ_IOINDEX(arg);
end);

#############################################################################
##
#F  PqProcessIndices() . . . . the list of active interactive ANUPQ processes
##
##  returns the list of (integer) indices of all active interactive  {\ANUPQ}
##  processes.
##
InstallGlobalFunction(PqProcessIndices, function()
  return Filtered( [1..Length(ANUPQData.io)], i -> IsBound( ANUPQData.io[i] ) );
end);

#############################################################################
##
#F  IsPqProcessAlive( <i> ) . .  checks an interactive ANUPQ process iostream
#F  IsPqProcessAlive()
##
##  return  `true'  if  the  {\GAP}  iostream  of  the  <i>th  (or   default)
##  interactive {\ANUPQ} process is alive (i.e. can still be written to),  or
##  `false', otherwise.
##
InstallGlobalFunction(IsPqProcessAlive, function(arg)
  return not IsEndOfStream( ANUPQData.io[ PqProcessIndex(arg) ].stream );
end);

#############################################################################
##
#V  PQ_MENUS . . . . . . . . . . . data describing the menus of the pq binary
##
##  a record whose fields are abbreviated names of  the  menus  of  the  `pq'
##  binary and whose values are themselves records with fields:
##
##    name
##        long name of menu;
##    depth
##        the number of times 0 must be passed to the `pq' binary for  it  to
##        exit;
##    prev
##        the menu one gets to from the current menu via option 0 (or `""' in
##        the case of the menu `SP';
##    nextopt
##        a record whose fields are the new menus of greater depth  that  can
##        be reached by an option of the current menu, and whose  values  are
##        the corresponding numbers of the options of the current menu needed
##        to get to the new menus.
##
InstallValue(PQ_MENUS, rec(
  SP  := rec( name  := "Standard Presentation Menu",
              depth := 1, prev  := "",   nextopt := rec( pQ := 7 ) ),
  pQ  := rec( name  := "(Main) p-Quotient Menu",
              depth := 2, prev  := "SP", nextopt := rec( pG  := 9, ApQ := 8 ) ),
  pG  := rec( name  := "(Main) p-Group Generation Menu",
              depth := 3, prev  := "pQ", nextopt := rec( ApG := 6 ) ),
  ApQ := rec( name  := "Advanced p-Quotient Menu",
              depth := 3, prev  := "pQ", nextopt := rec() ),
  ApG := rec( name  := "Advanced p-Group Gen'n Menu",
              depth := 4, prev  := "pG", nextopt := rec() )
  ) );

#############################################################################
##
#F  PQ_MENU( <datarec>, <newmenu> ) . . . . . . change/get menu of pq process
#F  PQ_MENU( <datarec> )
##
InstallGlobalFunction(PQ_MENU, function(arg)
local datarec, newmenu, nextmenu, tomenu, infolev;
  datarec := arg[1];
  if 2 = Length(arg) then
    newmenu := arg[2];
    if datarec.menu in ["SP", "pQ"] and newmenu in ["ApQ", "pG", "ApG"] then
      PQ_GRP_EXISTS_CHK( datarec ); #We try to avoid seg-faults!
    fi;
    while datarec.menu <> newmenu do
      if PQ_MENUS.(datarec.menu).depth >= PQ_MENUS.(newmenu).depth then
        datarec.menu := PQ_MENUS.(datarec.menu).prev;
        tomenu := PQ_MENUS.(datarec.menu).name;
        ToPQk(datarec, [ 0 ], [ "  #to ", tomenu]);
        infolev := 5;
      elif datarec.menu = "pQ" and newmenu = "ApQ" then
        datarec.menu := "ApQ";
        tomenu := PQ_MENUS.(datarec.menu).name;
        ToPQk(datarec, [ PQ_MENUS.pQ.nextopt.ApQ ], [ "  #to ", tomenu ]);
        infolev := 6;
      else
        nextmenu := RecNames( PQ_MENUS.(datarec.menu).nextopt )[1];
        tomenu := PQ_MENUS.(nextmenu).name;
        ToPQk(datarec, [ PQ_MENUS.(datarec.menu).nextopt.(nextmenu) ],
                       [ "  #to ", tomenu ]);
        datarec.menu := nextmenu;
        infolev := 6;
      fi;
      # menus are flushed at InfoANUPQ level 6, prompts at level 5
      if not IsBound(datarec.setupfile) then
        FLUSH_PQ_STREAM_UNTIL(datarec.stream, infolev, 5, PQ_READ_NEXT_LINE,
                              IS_PQ_PROMPT);
      fi;
    od;
  fi;
  return datarec.menu;
end);

#############################################################################
##
#F  IS_PQ_PROMPT( <line> ) . . . .  checks whether the line is a prompt of pq
##
##  returns `true' if the string  <line>  is  a  `pq'  prompt,  or  otherwise
##  returns `false'.
##
InstallGlobalFunction(IS_PQ_PROMPT,
  line -> IS_ALL_PQ_LINE(line) and ANUPQData.linetype = "prompt"
);

#############################################################################
##
#F  IS_ALL_PQ_LINE( <line> ) . checks whether line is a complete line from pq
##
##  returns `true' if the string <line> is a `pq' prompt or  a  request  from
##  `pq' to {\GAP} to compute stabilisers or simply ends  in  a  newline  and
##  sets `ANUPQData.linetype' to `"prompt"', `"request"'  or  `"hasnewline"',
##  accordingly; otherwise `ANUPQData.linetype' is  set  to  `"unknown"'  and
##  `false' is returned.
##
InstallGlobalFunction(IS_ALL_PQ_LINE, function( line )
local len;
  ANUPQData.linetype := "unknown";
  len := Length(line);
  if 0 < len then
    if line[len] = '\n' then
      if 4 < len  and line{[1 .. 3]} = "GAP" and line[len - 1] = '!' then
        ANUPQData.linetype := "request";
      elif 6 < len and line{[1 .. 6]} in ["Enter ", "Input "] then
        ANUPQData.linetype := "prompt";
      else
        ANUPQData.linetype := "hasnewline";
      fi;
    elif line = "Select option: " or
         1 < len and line{[len - 1 .. len]} = "? "  or
         8 < len and line{[len - 1 .. len]} = ": " and
                     line{[1 .. 6]} in ["Enter ", "Input ", "Add ne"] then
      ANUPQData.linetype := "prompt";
    fi;
  fi;
  return ANUPQData.linetype <> "unknown";
end);

#############################################################################
##
#F  PQ_READ_ALL_LINE( <iostream> ) .  read line from pq but poss. return fail
##
##  reads a complete line from <iostream> or return `fail'.
##
InstallGlobalFunction(PQ_READ_ALL_LINE,
  iostream -> ReadAllLine(iostream, false, IS_ALL_PQ_LINE)
);

#############################################################################
##
#F  PQ_READ_NEXT_LINE( <iostream> ) . read line from pq but never return fail
##
##  Essentially, like `PQ_READ_ALL_LINE' but we know there is a complete line
##  to be got, so we wait for it, before returning.
##
InstallGlobalFunction(PQ_READ_NEXT_LINE,
  iostream -> ReadAllLine(iostream, true, IS_ALL_PQ_LINE)
);

#############################################################################
##
#F  FLUSH_PQ_STREAM_UNTIL(<stream>,<infoLev>,<infoLevMy>,<readln>,<IsMyLine>)
##  . . .  . . . . . . . . . . . read lines from a stream until a wanted line
##
##  calls <readln> (which should be one of `ReadLine', `PQ_READ_NEXT_LINE' or
##  `PQ_READ_ALL_LINE') to read lines from a stream <stream> and `Info's each
##  line read at `InfoANUPQ' level <infoLev> until a line <line> is read  for
##  which `<IsMyLine>(<line>)' is `true'; <line> is `Info'-ed at  `InfoANUPQ'
##  level <infoLevMy> and returned. <IsMyLine>  should  be  a  boolean-valued
##  function that expects a string as its only argument,  and  <infoLev>  and
##  <infoLevMy> should be positive integers. An <infoLevMy> of 10 means  that
##  the  line  <line>  matched  by  `<IsMyLine>(<line>)'  should   never   be
##  `Info'-ed.
##
InstallGlobalFunction(FLUSH_PQ_STREAM_UNTIL,
function(stream, infoLev, infoLevMy, readln, IsMyLine)
local line;
  line := readln(stream);
  while not IsMyLine(line) do
    Info(InfoANUPQ, infoLev, Chomp(line));
    line := readln(stream);
  od;
  if line <> fail and infoLevMy < 10 then
    Info(InfoANUPQ, infoLevMy, Chomp(line));
  fi;
  return line;
end);

#############################################################################
##
#V  PQ_ERROR_EXIT_MESSAGES . . . error messages emitted by the pq before exit
##
##  A list of the error messages the `pq' emits just before exiting.
##
InstallValue(PQ_ERROR_EXIT_MESSAGES,
  [ "Evaluation in compute_degree may cause integer overflow",
    "A relation is too long -- increase the value of MAXWORD",
    "Ran out of space during computation" ]);

#############################################################################
##
#F  FILTER_PQ_STREAM_UNTIL_PROMPT( <datarec> )
##
##  reads `pq' output from `<datarec>.stream' until a `pq' prompt and `Info's
##  any lines that are prompts, blank lines, menu exits  or  start  with  the
##  strings in the list `<datarec>.filter' (if bound) at `InfoANUPQ' level 5;
##  all  other  lines  are  either  `Info'-ed  at  `InfoANUPQ'  level  3   if
##  `datarec.nonuser' is set, or, more usually, are `Info'-ed at  `InfoANUPQ'
##  level 2 if  they  are  computation  times  or  at  `InfoANUPQ'  level  1,
##  otherwise.
##
InstallGlobalFunction(FILTER_PQ_STREAM_UNTIL_PROMPT, function( datarec )
local match, filter, lowlev, ctimelev;
  filter := ["Exiting", "pq,", "Now enter",
             "Presentation listing images", "(use generators x1,x2"];
  if IsBound(datarec.match) then
    if datarec.match = true then
      match := ["Group:", "Group completed"];
    else
      match := [datarec.match];
    fi;
  fi;
  if IsBound(datarec.filter) then
    Append(filter, datarec.filter);
  fi;
  if ValueOption("nonuser") = true then
    lowlev := 3;
    ctimelev := 3;
  else
    ctimelev := 2;
    if not IsBound(datarec.OutputLevel) or datarec.OutputLevel = 0 then
      lowlev := 3;
    else
      lowlev := 1;
    fi;
  fi;
  repeat
    datarec.line := PQ_READ_NEXT_LINE(datarec.stream);
    if ANUPQData.linetype in ["prompt", "request"] then
      Info( InfoANUPQ, 5,        Chomp(datarec.line) );
      break;
    elif ForAny(["seconds", "Lused", "*** Final "],
                s -> PositionSublist(datarec.line, s) <> fail) then
      Info( InfoANUPQ, ctimelev, Chomp(datarec.line) );
    elif datarec.line = "\n" or
         ForAny( filter, s -> IsMatchingSublist(datarec.line, s) ) then
      Info( InfoANUPQ, 5,        Chomp(datarec.line) );
    elif PositionSublist(datarec.line, " saved on file") <> fail then
      Info( InfoANUPQ, ctimelev, Chomp(datarec.line) );
    elif ForAny( PQ_ERROR_EXIT_MESSAGES,
                 s -> IsMatchingSublist(datarec.line, s) ) then
      Info( InfoANUPQ + InfoWarning, 1, Chomp(datarec.line) );
      Error( "pq program terminated, with error condition:\n  ", datarec.line );
    else
      Info( InfoANUPQ, lowlev,   Chomp(datarec.line) );
    fi;
    if IsBound(match) then
      if ForAny( match, s -> IsMatchingSublist(datarec.line, s) ) then
        datarec.matchedline := datarec.line;
        datarec.complete := IsBound(datarec.complete) and datarec.complete or
                            IsMatchingSublist(datarec.line, "Group completed");
      fi;
    elif IsBound(datarec.matchlist) and
         ForAny( datarec.matchlist,
                 s -> PositionSublist(datarec.line, s) <> fail ) then
      Add(datarec.matchedlines, datarec.line);
    fi;
  until false;
end);

#############################################################################
##
#F  ToPQk( <datarec>, <cmd>, <comment> ) . . . . . . .  writes to a pq stream
##
##  writes  <cmd>  (and  <comment>,   in   setup   file   case)   to   stream
##  `<datarec>.stream' and `Info's <cmd> and <comment> at `InfoANUPQ' level 3
##  after a ```ToPQ> ''' prompt, and returns `true' if successful and  `fail'
##  otherwise. The ``k'' at the end of the  function  name  is  mnemonic  for
##  ``keyword'' (for ``keyword'' inputs to the `pq' binary one never wants to
##  flush output).
##
InstallGlobalFunction(ToPQk, function(datarec, cmd, comment)
local ok, line, i, j, closed, fragment, sepchars, words, filterones;

  if not IsOutputTextStream(datarec.stream) and
     IsEndOfStream(datarec.stream) then
    Error("sorry! Process stream has died!\n");
  fi;
  if cmd in ["gens", "rels"] then
    # these are done specially because of their potential to be enormously long
    if cmd = "gens" then
      line := "generators { ";
      sepchars := ", ";
    else
      line := "relators   { ";
      sepchars := "*^, ";
    fi;
    words := datarec.(cmd);
    filterones := cmd = "rels" and not IsBound(datarec.Relators) and
                  (IsFpGroup(datarec.group) or not IsPGroup(datarec.group));
    i := 1;
    while filterones and i <= Length(words) and IsOne(words[i]) do
      i := i + 1;
    od;
    if i <= Length(words) then
      Append(line, String(words[i]));
      i := i + 1;
    fi;
    ok := true;
    closed := false;
    repeat
      while filterones and i <= Length(words) and IsOne(words[i]) do
        i := i + 1;
      od;
      # i is the index of the next word to be added to line or > #words
      if i <= Length(words) then
        # if number of non-trivial words is 0 or 1 no comma is ever added
        Append(line, ", ");
        Append(line, String(words[i]));
        i := i + 1;
      else
        Append(line, " }");
        if cmd = "rels" then
          Append(line, ";");
        fi;
        closed := true; # not quite equivalent to: i > Length(words)
      fi;
      while ok and (Length(line) >= 69 or (closed and Length(line) > 0)) do
        if Length(line) >= 69 then
          # find a nice break if we can
          j := 68;
          while j > 4 and not line[j] in sepchars do j := j - 1; od;
          # no nice break
          if j = 4 then
            j := 69;
            while j < Length(line) and not line[j] in sepchars do
              j := j + 1;
            od;
          fi;
          fragment := line{[1 .. j]};
        else
          fragment := line;
          j := Length(line);
        fi;
        if j = Length(line) and closed then
          line := "";
        else
          line := Concatenation("  ", line{[j + 1 .. Length(line)]});
        fi;
        Info(InfoANUPQ, 4, "ToPQ> ", fragment);
        if IsBound( datarec.setupfile) then
          ok := WriteLine(datarec.stream, fragment);
        else
          ok := WriteLine(datarec.stream, fragment);
          if IsBound( ANUPQData.topqlogfile ) then
            WriteLine(ANUPQData.logstream, fragment);
          fi;
        fi;
      od;
    until closed or not ok;
  else
    # We add a null string in case <cmd> or <comment> is []
    # ... so that `Concatenation( List(., String) );' statements return strings
    Add(cmd, "");
    Add(comment, "");
    cmd     := Concatenation( List(cmd, String) );
    comment := Concatenation( List(comment, String) );
    Info(InfoANUPQ, 4, "ToPQ> ", cmd, comment);
    if IsBound( datarec.setupfile) then
      ok := WriteLine(datarec.stream, Concatenation(cmd, comment));
    else
      ok := WriteLine(datarec.stream, cmd);
      if IsBound( ANUPQData.topqlogfile ) then
        WriteLine(ANUPQData.logstream, Concatenation(cmd, comment));
      fi;
    fi;
  fi;
  if ok = fail then
    Error("write to stream failed\n");
  fi;
  return ok;
end);

#############################################################################
##
#F  ToPQ(<datarec>, <cmd>, <comment>) . .  write to pq (& for iostream flush)
##
##  calls `ToPQk' to write <cmd> (and  <comment>,  in  setup  file  case)  to
##  stream `<datarec>.stream' and `Info' <cmd> and <comment>  at  `InfoANUPQ'
##  level 3 after a ```ToPQ> ''' prompt, and then, if we are not just writing
##  a setup file (determined by  checking  whether  `<datarec>.setupfile'  is
##  bound), calls `FILTER_PQ_STREAM_UNTIL_PROMPT' to filter lines  to  `Info'
##  at the various `InfoANUPQ' levels. If we are not writing a setup file the
##  last line flushed is saved in `<datarec>.line'.
##
InstallGlobalFunction(ToPQ, function(datarec, cmd, comment)
  ToPQk(datarec, cmd, comment);
  if not IsBound( datarec.setupfile ) then
    FILTER_PQ_STREAM_UNTIL_PROMPT(datarec);

    while ANUPQData.linetype = "request" do
      HideGlobalVariables( "ANUPQglb", "F", "gens", "relativeOrders",
                           "ANUPQsize", "ANUPQagsize" );
      Read( Filename( ANUPQData.tmpdir, "GAP_input" ) );
      Read( Filename( ANUPQData.tmpdir, "GAP_rep" ) );
      UnhideGlobalVariables( "ANUPQglb", "F", "gens", "relativeOrders",
                             "ANUPQsize", "ANUPQagsize" );
      ToPQk( datarec, [ "pq, stabiliser is ready!" ], [] );
      FILTER_PQ_STREAM_UNTIL_PROMPT(datarec);
    od;
  fi;
end);

#############################################################################
##
#F  ToPQ_BOOL( <datarec>, <optval>, <comment> ) . . . .  pass a boolean to pq
##
##  converts a {\GAP} boolean  <optval>  to  a  C  boolean  and  appends  the
##  appropriate adjustment to the string <comment> before calling `ToPQ'  (we
##  assume that <optval> is boolean ... `VALUE_PQ_OPTION' should already have
##  checked that).
##
InstallGlobalFunction( ToPQ_BOOL, function( datarec, optval, comment )
  if optval = true then
    ToPQ( datarec, [ 1 ], [ "  #do ", comment ] );
  else
    ToPQ( datarec, [ 0 ], [ "  #do not ", comment ] );
  fi;
end);

#############################################################################
##
#F  PqRead( <i> )  . . .  primitive read of a single line from ANUPQ iostream
#F  PqRead()
##
##  read a complete line of  {\ANUPQ}  output,  from  the  <i>th  or  default
##  interactive {\ANUPQ} process, if there is output to be read  and  returns
##  `fail' otherwise. When successful, the  line  is  returned  as  a  string
##  complete with trailing newline, colon, or question-mark character. Please
##  note that it is possible to be ``too  quick''  (i.e.~the  return  can  be
##  `fail' purely because the output from {\ANUPQ} is not there yet), but  if
##  `PqRead' finds any output at all, it waits for a complete line.  `PqRead'
##  also writes the line read via `Info' at `InfoANUPQ' level 2.  It  doesn't
##  try to distinguish banner and menu output from other output of  the  `pq'
##  binary.
##
InstallGlobalFunction(PqRead, function(arg)
local line;

  line := PQ_READ_ALL_LINE( ANUPQData.io[ PqProcessIndex(arg) ].stream );
  Info(InfoANUPQ, 2, Chomp(line));
  return line;
end);

#############################################################################
##
#F  PqReadAll( <i> ) . . . . . read all current output from an ANUPQ iostream
#F  PqReadAll()
##
##  read and return as many *complete* lines of  {\ANUPQ}  output,  from  the
##  <i>th or default interactive {\ANUPQ} process, as there are to  be  read,
##  *at the time of the call*,  as  a  list  of  strings  with  any  trailing
##  newlines removed and returns the empty list otherwise.  `PqReadAll'  also
##  writes each line read via `Info' at `InfoANUPQ' level 2. It  doesn't  try
##  to distinguish banner and menu output  from  other  output  of  the  `pq'
##  binary. Whenever `PqReadAll' finds only a partial line, it waits for  the
##  complete line, thus increasing the probability that it has  captured  all
##  the output to be had from {\ANUPQ}.
##
InstallGlobalFunction(PqReadAll, function(arg)
local lines, stream, line;

  stream := ANUPQData.io[ PqProcessIndex(arg) ].stream;
  lines := [];
  line := PQ_READ_ALL_LINE(stream);
  while line <> fail do
    line := Chomp(line);
    Info(InfoANUPQ, 2, line);
    Add(lines, line);
    line := PQ_READ_ALL_LINE(stream);
  od;
  return lines;
end);

#############################################################################
##
#F  PqReadUntil( <i>, <IsMyLine> ) .  read from ANUPQ iostream until a cond'n
#F  PqReadUntil( <IsMyLine> )
#F  PqReadUntil( <i>, <IsMyLine>, <Modify> )
#F  PqReadUntil( <IsMyLine>, <Modify> )
##
##  read complete lines  of  {\ANUPQ}  output,  from  the  <i>th  or  default
##  interactive {\ANUPQ} process, ``chomps'' them (i.e.~removes any  trailing
##  newline character), emits them to `Info' at `InfoANUPQ' level 2  (without
##  trying to distinguish banner and menu output from  other  output  of  the
##  `pq' binary), and applies the function <Modify> (where <Modify>  is  just
##  the identity map/function for the first two forms)  until  a  ``chomped''
##  line  <line>  for  which  `<IsMyLine>(  <Modify>(<line>)  )'   is   true.
##  `PqReadUntil' returns the list of <Modify>-ed ``chomped'' lines read.
##
InstallGlobalFunction(PqReadUntil, function(arg)
local idx1stfn, stream, IsMyLine, Modify, lines, line;

  idx1stfn := First([1..Length(arg)], i -> IsFunction(arg[i]));
  if idx1stfn = fail then
    Error("expected at least one function argument\n");
  elif Length(arg) > idx1stfn + 1 then
    Error("expected 1 or 2 function arguments, not ",
          Length(arg) - idx1stfn + 1, "\n");
  elif idx1stfn > 2  then
    Error("expected 0 or 1 integer arguments, not ", idx1stfn - 1, "\n");
  else
    stream := ANUPQData.io[ ANUPQ_IOINDEX(arg{[1..idx1stfn - 1]}) ].stream;
    IsMyLine := arg[idx1stfn];
    if idx1stfn = Length(arg) then
      Modify := line -> line; # The identity function
    else
      Modify := arg[Length(arg)];
    fi;
    lines := [];
    repeat
      line := Chomp( PQ_READ_NEXT_LINE(stream) );
      Info(InfoANUPQ, 2, line);
      line := Modify(line);
      Add(lines, line);
    until IsMyLine(line);
    return lines;
  fi;
end);

#############################################################################
##
#F  PqWrite( <i>, <string> ) . . . . . . .  primitive write to ANUPQ iostream
#F  PqWrite( <string> )
##
##  write <string> to the <i>th  or  default  interactive  {\ANUPQ}  process;
##  <string> must be in exactly the form the {\ANUPQ} standalone expects. The
##  command is echoed via `Info' at `InfoANUPQ' level 3 (with a  ```ToPQ> '''
##  prompt); i.e.~do `SetInfoLevel(InfoANUPQ, 3);' to see what is transmitted
##  to the `pq' binary. `PqWrite' returns `true' if successful in writing  to
##  the stream of the interactive {\ANUPQ} process, and `fail' otherwise.
##
InstallGlobalFunction(PqWrite, function(arg)
local ioIndex, line;

  if Length(arg) in [1, 2] then
    ioIndex := ANUPQ_IOINDEX(arg{[1..Length(arg) - 1]});
    return ToPQk( ANUPQData.io[ioIndex], arg{[Length(arg)..Length(arg)]}, [] );
  else
    Error("expected 1 or 2 arguments ... not ", Length(arg), " arguments\n");
  fi;
end);

#############################################################################
##
#F  ANUPQ_ARG_CHK( <funcname>, <args> ) . . . . check args of int/non-int fns
##
##  checks the argument list <args> for a function that has both  interactive
##  and non-interactive versions, where <funcname> is the generic name of the
##  function. If <args> has length more than 1 then it contains  options  for
##  the function that have been passed in one of the {\GAP} 3-compatible ways
##  only available non-interactively. `ANUPQ_ARG_CHK' returns <datarec> which
##  is   either   `ANUPQData.ni'   in    the    non-interactive    case    or
##  `ANUPQData.io[<i>]' for some <i> in the interactive case,  after  setting
##  <datarec>.calltype' to one  of  `"interactive"',  `"non-interactive"'  or
##  `"GAP3compatible"'.
##
InstallGlobalFunction(ANUPQ_ARG_CHK, function(funcname, args)
local ioIndex, datarec, optrec, optnames;
  PQ_OTHER_OPTS_CHK( funcname, IsEmpty(args) or IsPosInt( args[1] ) );
  if IsEmpty(args) or IsPosInt( args[1] ) then
    datarec := ANUPQData.io[ CallFuncList( PqProcessIndex, args ) ];
    datarec.outfname := ANUPQData.outfile; # not always needed
    #datarec.calltype := "interactive";    # PqStart sets this
    if not IsBound(datarec.group) then
      Error( "huh! Interactive process has no group\n" );
    elif IsMatchingSublist(funcname, "PqDescendants") then
      if not IsPcGroup( datarec.group ) then
        Error( "group of process must be a pc group\n" );
      fi;
    else # Check for Prime, ClassBound if nec.
      PQ_OPTION_CHECK( funcname, datarec );
    fi;
  elif 1 = Length(args) then
    if not IsPcGroup( args[1] ) then
      if IsMatchingSublist(funcname, "PqDescendants") then
        Error( "first argument <args[1]> must be a pc group\n" );
      elif not IsFpGroup( args[1] ) then
        Error( "first argument <args[1]> must be a pc group or an fp group\n" );
      fi;
    fi;
    ANUPQData.ni := PQ_START( VALUE_PQ_OPTION( "PqWorkspace", 10000000 ),
                              VALUE_PQ_OPTION( "SetupFile" ) );
    datarec := ANUPQData.ni;
    datarec.group := args[1];
    datarec.calltype := "non-interactive";
    datarec.procId := 0;
    PQ_OPTION_CHECK( funcname, datarec ); # Check for Prime, ClassBound if nec.
    if IsBound( datarec.setupfile ) then
      datarec.outfname := "PQ_OUTPUT";
    else
      datarec.outfname := ANUPQData.outfile; # not always needed
    fi;
  else
    # GAP 3 way of passing options is supported in non-interactive use
    if funcname = "PqDescendantsTreeCoclassOne" then
      Error("GAP 3-compatible ways of passing options not supported");
    elif IsRecord(args[2]) then
      optrec := ShallowCopy(args[2]);
      optnames := Set( REC_NAMES(optrec) );
      SubtractSet( optnames, Set( ANUPQoptions.(funcname) ) );
      if not IsEmpty(optnames) then
        Error(ANUPQoptError( funcname, optnames ), "\n");
      fi;
    else
      optrec := ANUPQextractOptions(funcname, args{[2 .. Length(args)]});
    fi;
    PushOptions(optrec);
    PQ_FUNCTION.(funcname)( args{[1]} );
    PopOptions();
    datarec := ANUPQData.ni;
    datarec.calltype := "GAP3compatible";
    datarec.procId := 0;
  fi;
  return datarec;
end );

#############################################################################
##
#F  PQ_COMPLETE_NONINTERACTIVE_FUNC_CALL( <datarec> )
##
##  writes the final commands to the `pq' setup file so that the `pq'  binary
##  makes a clean exit, or just closes the stream to kill the `pq' process.
##
InstallGlobalFunction(PQ_COMPLETE_NONINTERACTIVE_FUNC_CALL, function(datarec)
  if IsBound( datarec.setupfile ) then
    PQ_MENU(datarec, "SP");
    ToPQk(datarec, [ 0 ], [ "  #exit program" ]);
  fi;
  CloseStream(datarec.stream);

  if IsBound( datarec.setupfile ) then
    Info(InfoANUPQ, 1, "Input file: '", datarec.setupfile, "' written.");
    Info(InfoANUPQ, 1, "Run `pq' with '", datarec.opts, "' flags.");
    Info(InfoANUPQ, 1, "The result will be saved in: '",
                       datarec.outfname, "'.");
  fi;
end );

#############################################################################
##
#F  ToPQLog([<filename>]) . . . . . . log or stop logging pq commands to file
##
##  With string argument <filename>,  `ToPQLog'  opens  the  file  with  name
##  <filename> for logging; all commands written to the `pq' binary (that are
##  `Info'-ed behind a ```ToPQ> ''' prompt at `InfoANUPQ' level 4)  are  then
##  also written to that  file  (but  without  prompts).  With  no  argument,
##  `ToPQLog' stops logging to whatever file was being logged to. If  a  file
##  was already being logged to, that file is closed and the file  with  name
##  <filename> is opened for logging.
##
InstallGlobalFunction(ToPQLog, function(arg)
  if not( IsEmpty(arg) or IsString( arg[1] ) ) then
    Error( "expected no arguments or one string argument\n" );
  fi;
  if IsBound(ANUPQData.topqlogfile) then
    CloseStream(ANUPQData.logstream);
    PQ_UNBIND(ANUPQData, ["topqlogfile", "logstream"]);
  elif IsEmpty(arg) then
    Info(InfoANUPQ + InfoWarning, 1, "No file currently being logged to.");
    return;
  fi;
  if not( IsEmpty(arg) ) and IsString(arg[1]) then
    ANUPQData.topqlogfile := arg[1];
    ANUPQData.logstream := OutputTextFile(ANUPQData.topqlogfile, false);
  fi;
end);

#E  anupqios.gi . . . . . . . . . . . . . . . . . . . . . . . . . . ends here


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

Class(CMScratchMainCodegen, ScratchMainCodegen, rec(
    genMain := (self, opts, sub, func1, func2, xy, params) >>
        func(TVoid, sub, Concatenation(xy, params), par_exec(
            ApplyFunc(call, Flat(Concatenation([func1],[false]))),
            ApplyFunc(call, Flat(Concatenation([func2], xy)))
        ))
));

Class(CMContextScratchCodegen, ScratchCodegen, rec(
    CPUCodegen := CPUCodegen,
    DMACodegen := DMACodegen,
    MainCodegen := CMScratchMainCodegen,
    Formula := meth(self, o, y, x, opts)
        local main, initfunc, cpufunc, dmafunc, auxfunc, prog, params, sub, initsub, memvars, scratchvars, v, io, loadvar, storevar, v, substrec, svars, initcode, datas, dvars, dv, tag;

        datas := Collect(o, FDataOfs);
        params := Set(Concatenation(Collect(o, param), Filtered(Collect(o, var), IsParallelLoopIndex)));
        sub := Cond(IsBound(opts.subName), opts.subName, "transform");

        initsub := Cond(IsBound(opts.subName), Concat("init_", opts.subName), "init");
        if IsBound(opts.generateInitFunc) and opts.generateInitFunc then
            initcode := chain(List(datas, x -> SReduce(x.var.init, opts)));
            initfunc := func(TVoid, initsub, params :: Set(Collect(initcode, param)), initcode);
        else
            initfunc := func(TVoid, initsub, params, chain());
        fi;

        dvars := List(datas, x->x.var);
        for dv in dvars do
            dv.t.qualifiers := [opts.romModifier];
        od;

        self.DMACodegen.loadbuffers := Set([]);
        self.DMACodegen.storebuffers := Set([]);
        self.DMACodegen.membuffers := Set([]);
        dmafunc := self.DMACodegen.Formula(o, y, x, opts);
        self.DMACodegen.membuffers := Set(self.DMACodegen.membuffers);
        SubtractSet(self.DMACodegen.membuffers, Set([x, y]));
        for v in self.DMACodegen.membuffers do v.t.qualifiers := [opts.memModifier]; od;

        self.CPUCodegen.loadbuffers := Set([]);
        self.CPUCodegen.storebuffers := Set([]);
        cpufunc := self.CPUCodegen.Formula(o, y, x, opts);

        memvars := Set(Concat(
            Filtered(Flat(List(Collect(cpufunc, decl), i->i.vars)), j->IsBound(j.t.qualifiers) and opts.memModifier in j.t.qualifiers),
            self.DMACodegen.membuffers));
        scratchvars := Set(Filtered(Flat(List(Collect(cpufunc, decl), i->i.vars)), j->IsBound(j.t.qualifiers) and opts.scratchModifier in j.t.qualifiers));

        tag := opts.tags[1];

        loadvar := var.fresh_t("S", TArray(opts.XType.t, Cond(tag.isRegCx, 2 * tag.size, tag.size)));
        loadvar.t.qualifiers := [opts.scratchModifier];
        storevar := var.fresh_t("S", TArray(opts.XType.t, Cond(tag.isRegCx, 2 * tag.size, tag.size)));
        storevar.t.qualifiers := [opts.scratchModifier];

        substrec := rec();
        for v in self.CPUCodegen.loadbuffers do substrec.(v.id) := loadvar; od;
        for v in self.CPUCodegen.storebuffers do substrec.(v.id) := storevar; od;
        for v in self.DMACodegen.loadbuffers do substrec.(v.id) := loadvar; od;
        for v in self.DMACodegen.storebuffers do substrec.(v.id) := storevar; od;
        cpufunc :=  SubstVars(cpufunc, substrec);
        dmafunc :=  SubstVars(dmafunc, substrec);

        cpufunc := SubstTopDown(cpufunc, @(1, decl), e->decl(Filtered(@(1).val.vars, i->not i in Concat(memvars, scratchvars, [loadvar, storevar])), @(1).val.cmd));
        dmafunc := SubstTopDown(dmafunc, @(1, decl), e->decl(Filtered(@(1).val.vars, i->not i in Concat(memvars, scratchvars, [loadvar, storevar])), @(1).val.cmd));

        [cpufunc, dmafunc, svars] := When(opts.swp,_double_buffer(cpufunc, dmafunc, [loadvar, storevar], opts),[cpufunc,dmafunc,[loadvar,storevar]]);

        x.t.qualifiers := Set(x.t.qualifiers);
        y.t.qualifiers := Set(y.t.qualifiers);

        main := self.MainCodegen.genMain(opts, sub, func(TVoid,"DMA_jump",[ false ],dma_jump(false)), dmafunc, [y,x], params);
        prog := program(
            decl(Concat(Set(memvars), svars, dvars), chain(
                initfunc,
                cpufunc,
		dmafunc,
                main
        )));

        prog.dimensions := o.dims();
        return prog;
    end
));


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

Class(phi, Exp);

Class(SSAMap, rec(
     __call__ := self >> WithBases(self, rec(
	         operations:=PrintOps,
		 top := 0 ,
		 stack := [ ],
		 alternatives := [ ],
		 phi_map := tab()
	     )).enter(),
     print := self >> Print("SSAMap(", self.top, ")"),

     insert := meth(self, oldvar, newvar)
         self.stack[self.top].(oldvar.id) := newvar;
     end,

     _lookup := meth(self, v, i)
	 while i >= 1 do
	     if IsBound(self.stack[i].(v)) then return self.stack[i].(v); fi;
	     i := i-1;
	 od;
	 return false;
     end,

     lookup := (self, v) >> self._lookup(v.id, self.top),

     enter := meth(self)
         local newtab;
	 newtab := tab();
         Add(self.stack, newtab);
	 Add(self.alternatives, [newtab]);
	 self.top := self.top + 1;
	 return self;
     end,

     next := meth(self)
         local newtab;
	 newtab := tab();
         Add(self.alternatives[self.top], newtab);
	 self.stack[self.top] := newtab;
	 return self;
     end,

     join := meth(self)
         local i, alt, tt, id, phi_expr, phi_var, outer_tab, type, alts, tmp;

     #Create a list of all the changed variables in merging scope
#      changed := tab();

#      for alt in self.alternatives[self.top] do
#          for id in Filtered(NSFields(alt),x->x<>"__doc__") do
#              if not(IsBound(changed.(id))) then
#                  changed.(id) := (id);
#              fi;
#          od;
#      od;

     #add alternatives for when variable appear in two different branches
	 tt := tab();
	 for alt in self.alternatives[self.top] do
	     for id in Filtered(NSFields(alt),x->x<>"__doc__") do
	         if IsBound(tt.(id)) then Add(tt.(id).args, alt.(id));
		 else
		     tt.(id) := phi(alt.(id));
		     tt.(id).t := (alt.(id).t);
		 fi;
	     od;
	 od;

     #if a variable does not appear in all branches, maybe it should be set do its
     #father version
     alts := Length(self.alternatives[self.top]);
     for id in Filtered(NSFields(tt),x->x<>"__doc__") do
         if (Length(tt.(id).args)<>alts)  then
             tmp := self._lookup((id), self.top-1);
             if tmp<>false then
                 Add(tt.(id).args, tmp);
             fi;
         fi;
     od;

	 outer_tab := self.stack[self.top-1];
	 self.phi_map := tab(); # phiXX -> phi(a, b, c)
	 for id in Filtered(NSFields(tt),x->x<>"__doc__") do
	     phi_expr := tt.(id);
	     if Length(phi_expr.args) > 1 then # non-trivial phi function
		 phi_var := var.fresh_t("phi", phi_expr.args[1].t);
                 phi_expr.t := phi_var.t;
		 self.phi_map.(phi_var.id) := phi_expr;
		 outer_tab.(id) := phi_var;
		 # we make the following assumptions here
		 # 1. variables inside phi(...) are live outside the if
		 #    (violated when result of phi(...) is never used, eg. phi itself is dead,
		 #     this is safe, but will not remove dead code)
		 # 2. variables not inside phi(...) are dead outside the if
		 #    (violated in obscure cases like if(x) {t1=1} else {}; if(x) {y=t1}else{y=1},
		 #     we can prove that these cases do not arise in linear transforms)
		 for i in phi_expr.args do
		     i.live_out := true;
		 od;
	     else  # trivial phi(x)->x
             outer_tab.(id) := phi_expr.args[1];
	     fi;
	 od;
	 return self;
     end,

     leave := meth(self)
         local phi_map;
         phi_map := self.join();
         RemoveLast(self.stack, 1);
         RemoveLast(self.alternatives, 1);
	 self.top := self.top - 1;
	 return phi_map;
     end,

     subst_vars := (self, exp) >>
         FoldR(self.stack, (ex, t) -> SubstVars(ex, t), exp),

     phi_map_assigns := self >> List(Filtered(NSFields(self.phi_map),x->x<>"__doc__"),
	 f -> assign(var(f), self.phi_map.(f)))

));

Declare(_SSA);

SSA := function(code)
    local ssamap, icode;
    ssamap := SSAMap();
    icode := FlattenCode(_SSA(code, ssamap));
    return icode;
end;

_SSA := function(code, ssamap)
    local i, c, newvar, newc, phi_map;

    if ObjId(code) = IF then
	code.cond := ssamap.subst_vars(code.cond);
	ssamap.enter();	code.then_cmd := _SSA(code.then_cmd, ssamap);
	ssamap.next();	code.else_cmd := _SSA(code.else_cmd, ssamap);
	ssamap.leave();
	return chain(code, ssamap.phi_map_assigns());

    elif ObjId(code) = assign then
	return _SSA(chain(code), ssamap);

    elif IsExpCommand(code) then
	return ssamap.subst_vars(code);

    elif ObjId(code) = chain then
	for i in [1..Length(code.cmds)] do
	    c := code.cmds[i];
	    if ObjId(c) = assign then
		    c.exp := ssamap.subst_vars(c.exp);
		    if IsVar(c.loc) then
			newvar := var.fresh_t("t", c.loc.t);
			ssamap.insert(c.loc, newvar);
			c.loc := newvar;
		    else
			c.loc := ssamap.subst_vars(c.loc);
		    fi;
	    else
		code.cmds[i] := _SSA(c, ssamap);
	    fi;
	od;
	return code;
    elif IsRec(code) and IsBound(code.sideeffect) and code.sideeffect then
        return map_children_safe(code, x->ssamap.subst_vars(x));
    elif IsCommand(code) then
        return code;
    else
        Error("Can't handle ", ObjId(code));
    fi;
 end;

StronglyConnectedComponents := function(ll)
    local result, l, bool, r;

    result := [];
    for l in ll do
        bool := true;
        for r in [1..Length(result)] do
            if (bool=true) then
                if Length(Intersection(l, result[r]))<>0 then
                    result[r]:=Concatenation(result[r], l);
                    bool:=false;
                fi;
            fi;
        od;
        if (bool=true) then
            Add(result, l);
        fi;
    od;

    return result;
end;

EliminatePhiSSA := function(code)
    local elim_map, phis, p, v, master_var;
    elim_map := tab();
    phis := Collect(code, phi);

    phis := StronglyConnectedComponents(List(phis, x->x.args));

    for p in phis do
        master_var := p[1];
        for v in Drop(p, 1) do
	        elim_map.(v.id) := master_var;
    	od;
    od;
    code := SubstVars(code, elim_map);
    code := SubstBottomUp(code, phi, e->e.args[1]);
    return code;
end;


text:
--
-- An instance of this template must have a %Export section and the export_terminals option
--
-- Macros that may be redefined in an instance of this template
--
--     %eof_token
--     %additional_interfaces
--     %super_stream_class -- subclass com.ibm.lpg.Utf8LpgLexStream for getKind
--     %prs_stream_class -- use /.PrsStream./ if not subclassing
--
-- B E G I N N I N G   O F   T E M P L A T E   LexerTemplateD
--
%Options programming_Language=typescript,margin=4
%Options table
%options action-block=("*.ts", "/.", "./")
%options ParseTable=ParseTable
%Options prefix=Char_

--
-- This template requires that the name of the EOF token be set
-- to EOF and that the prefix be "Char_" to be consistent with
-- KeywordTemplateD.
--
%Eof
    EOF
%End

--
-- This template also requires that the name of the parser EOF
-- Token to be exported be set to EOF_TOKEN
--
%Export
    EOF_TOKEN
%End

%Define
    --
    -- Macros that are be needed in an instance of this template
    --
    $eof_token /.%_EOF_TOKEN./

    $additional_interfaces /../
    $super_stream_class /.%file_prefix%Utf8LpgLexStream./
    $prs_stream_class /.IPrsStream./
    $super_class /.Object./

    --
    -- Macros useful for specifying actions
    --
    $Header
    /.
                //
                // Rule %rule_number:  %rule_text
                //
                ./

    $DefaultAction
    /.%Header%case %rule_number: { ./

    $BeginAction /.%DefaultAction./

    $EndAction
    /.          break;
                }./

    $BeginJava
    /.%BeginAction
                %symbol_declarations./

    $EndJava /.%EndAction./

    $NoAction
    /.%Header%case %rule_number:
                    break; ./

    $BeginActions
    /.
        public void ruleAction( ruleNumber : number )
        {
            switch(ruleNumber)
            {./

    $SplitActions
    /.
	            default:
	                ruleAction%rule_number(ruleNumber);
	                break;
	        }
	        return;
	    }

	    public void ruleAction%rule_number(ruleNumber : number )
	    {
	        switch (ruleNumber)
	        {./

    $EndActions
    /.
                default:
                    break;
            }
            return;
        }./
%End

%Globals
    /.
    import { RuleAction, ParseTable, LexParser, ILexStream, IPrsStream, Monitor, LpgLexStream } from "lpg2ts";
    ./
%End

%Headers
    /.
    export class %action_type extends %super_class implements RuleAction%additional_interfaces
    {
        private %super_stream_class utf8LexStream;

        private static ParseTable prs = new %prs_type();
        public ParseTable getParseTable() { return prs; }

        private LexParser lexParser = new LexParser();
        public LexParser getParser() { return lexParser; }

        public number getToken(i : number) { return lexParser.getToken(i); }
        public number getRhsFirstTokenIndex(i : number) { return lexParser.getFirstToken(i); }
        public number getRhsLastTokenIndex(i : number) { return lexParser.getLastToken(i); }

        public number getLeftSpan() { return lexParser.getToken(1); }
        public number getRightSpan() { return lexParser.getLastToken(); }

        public void resetKeywordLexer()
        {
            if (!this.kwLexer)
                  this.kwLexer = new %kw_lexer_class(utf8LexStream.getInputBytes(), %_IDENTIFIER);
            else this.kwLexer.setInputBytes(utf8LexStream.getInputBytes());
        }

        public void reset(filename : string, number tab)
        {
            utf8LexStream = new %super_stream_class(filename, tab);
            lexParser.reset((ILexStream) utf8LexStream, prs, (RuleAction) this);
            resetKeywordLexer();
        }

        public void reset(byte[] input_bytes, filename : string)
        {
            reset(input_bytes, filename, 1);
        }

        public void reset(byte[] input_bytes, filename : string, number tab)
        {
            utf8LexStream = new %super_stream_class(input_bytes, filename, tab);
            lexParser.reset((ILexStream) utf8LexStream, prs, (RuleAction) this);
            resetKeywordLexer();
        }

        public %action_type(filename : string, number tab)
        {
            reset(filename, tab);
        }

        public %action_type(byte[] input_bytes, filename : string, number tab)
        {
            reset(input_bytes, filename, tab);
        }

        public %action_type(byte[] input_bytes, filename : string)
        {
            reset(input_bytes, filename, 1);
        }

        public %action_type() {}

        public ILexStream getILexStream() { return utf8LexStream; }

        /**
         * @deprecated replaced by {@link #getILexStream()}
         */
        public ILexStream getLexStream() { return utf8LexStream; }

        private void initializeLexer(%prs_stream_class this.prsStream, number start_offset, number end_offset)
        {
            if (utf8LexStream.getInputBytes() == null)
                throw new ReferenceError("LexStream was not initialized");
            utf8LexStream.setPrsStream(this.prsStream);
            this.prsStream.makeToken(start_offset, end_offset, 0); // Token list must start with a bad token
        }

        public void lexer(%prs_stream_class this.prsStream)
        {
            lexer(null, this.prsStream);
        }

        public void lexer(Monitor monitor, %prs_stream_class this.prsStream)
        {
            if (utf8LexStream.getInputBytes() == null)
                throw new ReferenceError("Utf8LexStream was not initialized");

            utf8LexStream.setPrsStream(this.prsStream);

            this.prsStream.makeToken(0, 0, 0); // Token list must start with a bad token

            lexParser.parseCharacters(monitor);  // Lex the input characters

            i : number = utf8LexStream.getStreamIndex();
            this.prsStream.makeToken(i, i, %eof_token); // and end with the end of file token
            this.prsStream.setStreamLength(this.prsStream.getSize());

            return;
        }
    ./
%End

%Rules
    /.%BeginActions./
%End

%Trailers
    /.
        %EndActions
    }
    ./
%End

--
-- E N D   O F   T E M P L A T E
--