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

TrPRDFT1 := arg -> let(x:=ApplyFunc(PRDFT1, arg).transpose(), Chain(PrintLine(x, " ", x.dimensions), x));
TrPRDFT2 := arg -> ApplyFunc(PRDFT2, arg).transpose();
TrPRDFT3 := arg -> ApplyFunc(PRDFT3, arg).transpose();
TrPRDFT4 := arg -> ApplyFunc(PRDFT4, arg).transpose();

DFTSymmetries := rec(    #n-even, #n-odd
    DFT  := [
	[Real ,                 [CE0_R2Cpx(W_RFT, 1), [PRDFT1,   PRDFT1], 1]],
	[CE0_R2Cpx(W_URFTT, 1), [Real,                [TrPRDFT1, TrPRDFT1], 1]],
	[Even0, [Even0, [DCT1, DCT5], 1]],
	[Odd00, [Odd00, [DST1, DST5], E(4)]]  # *j
    ],

    DFT2 := [
	[Real ,                 [CO0_R2Cpx(W_RFT, -E(4)), [PRDFT2,   PRDFT2], 1]],
	[CE1_R2Cpx(W_URFTT, 1), [Real,                    [TrPRDFT3, TrPRDFT3], 1]],
	[Even1, [Odd0,   [DCT2, DCT6], 1]],
	[Odd1 , [Even00, [DST2, DST6], E(4)]]
    ],

    DFT3 := [
	[Real ,                    [CE1_R2Cpx(W_RFT, 1), [PRDFT3,   PRDFT3], 1]],
	[CO1_R2Cpx(W_URFTT, E(4)), [Real,                [TrPRDFT2, TrPRDFT2], 1]],
	[Real , [CE1_R2Cpx(W_RFT, 1), [PRDFT3, PRDFT3], 1]],
	[Odd0 , [Even1, [DCT3, DCT7], 1]],
	[Even00,[Odd1,  [DST3, DST7], E(4)]]
    ],

    DFT4 := [
	[Real ,                    [CO1_R2Cpx(W_RFT, -E(4)), [PRDFT4,   PRDFT4], 1]],
	[CO1_R2Cpx(W_URFTT, E(4)), [Real,                    [TrPRDFT4, TrPRDFT4], 1]],
	[Odd1 , [Odd1,  [DCT4, DCT8], 1]],
	[Even1, [Even1, [DST4, DST8], E(4)]]
    ]
);

_dftSym := function(nt, sym)
    local pos, entry, dftsym;
    if not IsBound(DFTSymmetries.(nt.name))
	then return false; fi;

    dftsym := DFTSymmetries.(nt.name);
    pos := PositionProperty(dftsym, x->x[1]=sym);
    if pos = false then
	return false; fi;

    entry := dftsym[pos][2];
    if IsEvenInt(Rows(nt)) then return [entry[1], entry[2][1], entry[3]];
    else                        return [entry[1], entry[2][2], entry[3]];
    fi;
end;

# let Q := JDFTp_Reconstruct(n), p = 1..4 (DFT type)
# Q is defined to be the matrix that satisfies
#    Q * DFTp(n) = DFTp(n) * J(n)
#    Q = J(n) ^ (DFTp(n)^-1)
#
#spiral> PrintMat( MatSPL(J(8))^(MatSPL(DFT(8))^-1) );

Class(JRecTwid1, DiagFunc, rec(
    def := n -> rec(size:=n),
    range := self >> TComplex,
    lambda := self >> let(n:=self.params[1], i:=Ind(n),
	Lambda(i, omega(n, -i)))
));

Class(JRecTwid3, DiagFunc, rec(
    def := n -> rec(size:=n),
    range := self >> TComplex,
    lambda := self >> let(n:=self.params[1], i:=Ind(n),
	Lambda(i, omega(2*n, (-2*i + n - 1))))
));

mk_InvDFT_Reconstruct := dft -> (let(n:=dft.params[1], Cond(
    ObjId(dft) = DFT1, DirectSum(I(1), J(n-1)),
    ObjId(dft) = DFT2, DirectSum(I(1), -J(n-1)),
    ObjId(dft) = DFT3, J(n),
    ObjId(dft) = DFT4, -J(n),
    Error("<dft> must be a non-terminal DFT1, DFT2, DFT3, or DFT4"))));

mk_JInvDFT_Reconstruct := dft -> (let(n:=dft.params[1], Cond(
    ObjId(dft) = DFT1, Diag(JRecTwid1(n)),
    ObjId(dft) = DFT2, I(n),
    ObjId(dft) = DFT3, Diag(JRecTwid3(n)),
    ObjId(dft) = DFT4, -I(n),
    Error("<dft> must be a non-terminal DFT1, DFT2, DFT3, or DFT4"))));

mk_JDFT_Reconstruct := dft -> (let(n:=dft.params[1], Cond(
    ObjId(dft) = DFT1, Diag(JRecTwid1(n)) * mk_InvDFT_Reconstruct(dft),
    ObjId(dft) = DFT2, mk_InvDFT_Reconstruct(dft),
    ObjId(dft) = DFT3, Diag(JRecTwid3(n)) * mk_InvDFT_Reconstruct(dft),
    ObjId(dft) = DFT4, -mk_InvDFT_Reconstruct(dft),
    Error("<dft> must be a non-terminal DFT1, DFT2, DFT3, or DFT4"))));

verify_JInvDFT_Reconstruct := dft -> let(T:=ObjId(dft), n:=dft.params[1], k:=dft.params[2],
    MatSPL(dft * J(n)) - MatSPL(mk_JInvDFT_Reconstruct(dft) * T(n, -k))
);
verify_JDFT_Reconstruct := dft -> let(T:=ObjId(dft), n:=dft.params[1], k:=dft.params[2],
    MatSPL(dft * J(n)) - MatSPL(mk_JDFT_Reconstruct(dft) * dft)
);

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

     DFT1_upgrade := ARule(Compose, [@(1,DFT), [@(2,GathExtend), @, @(3).cond(e->e in [Even1, Odd1])]],
	 e -> let(n := Cols(@(1).val), [ Diag(Twid(n,n,-1,1/2,0,0)), DFT2(n), @(2).val])),

     DFT3_upgrade := ARule(Compose, [@(1,DFT3), [@(2,GathExtend), @, @(3).cond(e->e in [Even1, Odd1])]],
	 e -> let(n := Cols(@(1).val), [ Diag(Twid(n,n,-1,1/2,1/2,0)), DFT4(n), @(2).val])),

     DFT1_upgradeCE := ARule(Compose, [@(1,DFT), [@(2,GathExtend), @, @(3, [CE1_R2Cpx, CO1_R2Cpx])]],
	 e -> let(n := Cols(@(1).val), [ Diag(Twid(n,n,-1,1/2,0,0)), DFT2(n), @(2).val])),

     DFT3_upgradeCE := ARule(Compose, [@(1,DFT3), [@(2,GathExtend), @, @(3, [CE1_R2Cpx, CO1_R2Cpx])]],
	 e -> let(n := Cols(@(1).val), [ Diag(Twid(n,n,-1,1/2,1/2,0)), DFT4(n), @(2).val])),

     DFT_EO_Symmetry := ARule(Compose,
	 [ @(1,[DFT,DFT2,DFT3,DFT4]),
	   [ @(2, GathExtend), @, @(3).cond(e->_dftSym(@(1).val, e)<>false)] ],
	 e ->
	     let(sym := _dftSym(@(1).val, @(3).val),
		 out := sym[1],
		 transf := sym[2],
		 n := When(@(3).val.tsize = "cols", Cols(@(2).val), Rows(@(2).val)),
		 scale := sym[3],
		 [PushL(GathExtendU(Rows(@(1).val), out)), scale*transf(n)])),

     L_Real := ARule(Compose, [[L, @(1), @(2)], [@(3,GathExtend), @, Real]],
	 e -> let(N:=@(1).val, k:=@(2).val,
	     [ GathExtend(N, Real), L(N, k) ])),

     L_Conj0_R2Cpx := ARule(Compose, [[L, @(1), @(2)], [@(3,GathExtend), @, @(4, [CE0_R2Cpx, CO0_R2Cpx])]],
	 e -> let(N:=@(1).val, k:=@(2).val, m:=N/k, nn:=Cols(@(3).val)/2,
	          sym  := @(4).val,
	          sym1 := When(ObjId(sym)=CE0_R2Cpx, CE1_R2Cpx(sym.w, sym.j), CO1_R2Cpx(sym.w, sym.j)),
		  diag := When(ObjId(sym)=CE0_R2Cpx, Diag(BHD(m, 1, -1)), Diag(BHD(m, -1, 1))),
	     [ DirectSum(GathExtend(m, sym), GathExtend(N-m, sym1)),
	       DirectSum(I(Int(m+1+_even(m))), Tensor(I(Int((k-1)/2)), diag), I(_even(k)*(m+_odd(m)))),
	       RC(VStack(Gath(Refl(nn, N, Int((m+2)/2), L(N,k))),
		         Gath(Refl(nn, N, nn-Int((m+2)/2), fCompose(L(N,k), fAdd(N, N-m, m)))))) ])),

     L_Conj1_R2Cpx := ARule(Compose, [[L, @(1), @(2)], [@(3,GathExtend), @, @(4, [CE1_R2Cpx, CO1_R2Cpx])]],
	 e -> let(N:=@(1).val, k:=@(2).val, m:=N/k, nn:=Cols(@(3).val)/2,
	          sym  := @(4).val,
		  diag := When(sym=CE1_R2Cpx, Diag(BHD(m, 1, -1)), Diag(BHD(m, -1, 1))),
	     [ GathExtend(N, sym),
	       DirectSum(Tensor(I(Int(k/2)), diag), When(IsOddInt(k), I(m+_odd(m)), [])),
	       RC(Gath(Refl(nn, N-1, nn, L(N,k)))) ])),

     L_Even0 := ARule(Compose, [[L, @(1), @(2)], [@(3,GathExtend), @, Even0]],
	 e -> let(N:=@(1).val, k:=@(2).val, m:=N/k, nn:=Cols(@(3).val),
	     [ DirectSum(GathExtend(m, Even0), GathExtend(N-m, Even1)),
	       VStack(Gath(Refl(nn, N, Int((m+2)/2), L(N,k))),
		      Gath(Refl(nn, N, nn-Int((m+2)/2), fCompose(L(N,k), fAdd(N, N-m, m))))) ])),

     L_Even00 := ARule(Compose, [[L, @(1), @(2)], [@(3,GathExtend), @, Even00]],
	 e -> let(N:=@(1).val, k:=@(2).val, m:=N/k, nn:=Cols(@(3).val),
	     [ DirectSum(GathExtend(m, Even00), GathExtend(N-m, Even1)),
	       VStack(Gath(Refl(nn, N-2, Int(m/2), fCompose(fAdd(N,N,-1), L(N,k), fAdd(N, N-1, 1)))),
		      Gath(Refl(nn, N-2, nn-Int(m/2), fCompose(fAdd(N,N,-1), L(N,k), fAdd(N, N-m, m))))) ])),

     L_Odd0 := ARule(Compose, [[L, @(1), @(2)], [@(3,GathExtend), @, Odd0]],
	 e -> let(N:=@(1).val, k:=@(2).val, m:=N/k, nn:=Cols(@(3).val),
	     [ DirectSum(GathExtend(m, Odd0),  GathExtend(m*(k-1), Odd1)),
	       DirectSum(I(Int((m+1)/2)), Tensor(I(Int((k-1)/2)), Diag(BHN(m))), When(IsEvenInt(k), I(Int(m/2)), [])),
	       VStack(Gath(Refl(nn, N, Int((m+1)/2), L(N,k))),
		      Gath(Refl(nn, N, nn-Int((m+1)/2), fCompose(L(N,k), fAdd(N, N-m, m))))) ])),

     L_Odd00 := ARule(Compose, [[L, @(1), @(2)], [@(3,GathExtend), @, Odd00]],
	 e -> let(N:=@(1).val, k:=@(2).val, m:=N/k, nn:=Cols(@(3).val),
	     [ DirectSum(GathExtend(m, Odd00), GathExtend(m*(k-1), Odd1)),
	       DirectSum(I(Int((m-1)/2)), Tensor(I(Int((k-1)/2)), Diag(BHN(m))), When(IsEvenInt(k), I(Int(m/2)), [])),
	       VStack(Gath(Refl(nn, N-2, Int((m-1)/2), fCompose(fAdd(N,N,-1), L(N,k), fAdd(N, N-1, 1)))),
		      Gath(Refl(nn, N-2, nn-Int((m-1)/2), fCompose(fAdd(N,N,-1), L(N,k), fAdd(N, N-m, m))))) ])),

     L_Even1 := ARule(Compose, [[L, @(1), @(2)], [@(3,GathExtend), @, Even1]],
	 e -> let(N:=@(1).val, k:=@(2).val, m:=N/k, nn:=Cols(@(3).val),
	     [ GathExtend(N, Even1),
	       Gath(Refl(nn, N-1, nn, L(N,k))) ])),

     L_Odd1 := ARule(Compose, [[L, @(1), @(2)], [@(3,GathExtend), @, Odd1]],
	 e -> let(N:=@(1).val, k:=@(2).val, m:=N/k, nn:=Cols(@(3).val),
	     [ GathExtend(N, Odd1),
	       DirectSum(Tensor(I(Int(k/2)), Diag(BHN(m))), When(IsOddInt(k), I(Int(m/2)), [])),
	       Gath(Refl(nn, N-1, nn, L(N,k))) ])),

     CRT_Even1 := ARule(Compose, [@(1, CRT), [@(2,GathExtend), @, Even1]],
	 e -> let(N:=Rows(@(2).val),  nn:=Cols(@(2).val),
	     [ GathExtend(N, Even1),
	       Gath(Refl(nn, N-1, nn, @(1).val)) ])),

     Tensor_Dirsum_split := ARule(Compose,
	 [[Tensor, [I, @(1)], @(2)], [DirectSum, @(3).cond(e->Rows(e)=Cols(@(2).val)), @(4)]],
	 e -> [ DirectSum(@(2).val * @(3).val, Tensor(I(@(1).val-1), @(2).val) * @(4).val) ]),

     IxDFT_EO1 := ARule(Compose,
	 [[Tensor, [I, @(1)], [@(2,[DFT,DFT2,DFT3,DFT4]), @, 1, ...]], [@(3,GathExtend), @, @(4).cond(e->e in [Even1, Odd1])]],
	 e -> let(k:=@(1).val, n:=Cols(@(2).val), transf:=@(2).val, sign := When(@(4).val=Odd1, -1, 1),
	        When(k mod 2 = 0,
	         [ # Below is actually a scaled type 1 extension! exploit this for final rewriting phase
                   PushL(DirectSum(I(n*k/2), Tensor(I(k/2), mk_JDFT_Reconstruct(transf)*J(n))) *
                         @(3).val) *
		   Tensor(I(k/2), transf) ],

                   # again, a scaled type 1 extension!
		 [ PushL(VStack(I(n*(k+1)/2),
			        Tensor(J((k-1)/2), sign * mk_JDFT_Reconstruct(transf))
				   * Gath(fAdd(n*(k+1)/2, n*(k-1)/2, 0)))),
		   DirectSum(Tensor(I((k-1)/2), 1/2*transf),
			     transf * GathExtend(n, @(4).val))]))),

 #     IxDFT_EO1 := ARule(Compose,
# 	 [[Tensor, [I, @(1)], [@(2,[DFT,DFT2,DFT3,DFT4]), @, 1, ...]], [GathExtend, @(3), @(4).cond(e->e in [Even1, Odd1])]],
# 	 e -> let(k:=@(1).val, n:=Cols(@(2).val), transf:=@(2).val, sign := When(@(4).val=Odd1, -1, 1),
# 	        When(k mod 2 = 0,
# 	         [ # Below is actually a scaled type 1 extension! exploit this for final rewriting phase
#                    PushL(VStack(I(n*k/2), Tensor(J(k/2), sign * mk_JDFT_Reconstruct(transf)))),
# 		   Tensor(I(k/2), 1/2*transf) ],

#                    # again, a scaled type 1 extension!
# 		 [ PushL(VStack(I(n*(k+1)/2),
# 			        Tensor(J((k-1)/2), sign * mk_JDFT_Reconstruct(transf))
# 				   * Gath(fAdd(n*(k+1)/2, n*(k-1)/2, 0)))),
# 		   DirectSum(Tensor(I((k-1)/2), 1/2*transf),
# 			     transf * GathExtend(n, @(4).val))]))),

     IxDFT_CEO1_RFTT_DHT := ARule(Compose,
	 [[Tensor, [I, @(1)], @(2,[DFT,DFT2,DFT3,DFT4])],
	  [GathExtend, @(3), @(4,[CE1_R2Cpx, CO1_R2Cpx],e->e.w in [W_RFTT, W_URFTT, W_DHT, W_DHT])]],

	 e -> let(k:=@(1).val, n:=Cols(@(2).val),
	          transf := @(2).val, sign := When(ObjId(@(4).val)=CO1_R2Cpx, -1, 1),
	          newtransf := RC(ObjId(transf)(transf.params[1], -transf.params[2])),
		  winv := @(4).val.w^-1,
		  pack1 := Mat([winv[1]]),
		  pack2 := Mat([winv[2]]),

	        When(k mod 2 = 0,
	         [ PushL(VStack(Tensor(I(n*k/2), pack1),
			        Tensor(J(k/2), sign * Tensor(I(n), pack2)*
				               RC(mk_JInvDFT_Reconstruct(transf))))),
		   Tensor(I(k/2), newtransf) ],

		 [ PushL(VStack(DirectSum(Tensor(I(n*(k-1)/2), pack1), I(n)),
			        Tensor(J((k-1)/2), sign * Tensor(I(n), pack2) *
				                   RC(mk_JInvDFT_Reconstruct(transf))) * Gath(fAdd(n*k, n*(k-1), 0)))),
		   DirectSum(Tensor(I((k-1)/2), newtransf),
			     transf * GathExtend(n, @(4).val))]))),

     IxDFT_Real := ARule(Compose,
	 [[Tensor, [I, @(1)], [@(2,[DFT,DFT2,DFT3,DFT4]), @, 1, ...]], [GathExtend, @(3), Real]],
	 e -> [ Tensor(I(@(1).val), @(2).val * GathExtend(Cols(@(2).val), Real)) ]),

     DFTxI_Ext := ARule(Compose,
	 [[@(1,Tensor), [@(2,[DFT,DFT2,DFT3,DFT4]), @, 1, ...], I], @(4, GathExtend)],
	 e -> [ @(1).val.parallelForm() * @(4).val]),

     Tensor_J1 := Rule([Tensor, [J, 1], @(1)], e -> @(1).val),
     Tensor_I1 := Rule([Tensor, [I, 1], @(1)], e -> @(1).val),
     Tensor_I0 := Rule([Tensor, ..., [I, 0], ...], e -> I(0)),
     Compose_I := ARule(Compose, [I], e -> []),
     DirectSum_I0 := ARule(DirectSum, [[I, 0]], e -> []),
     DirectSum_Assoc := ARule(DirectSum, [@(1,DirectSum)], e -> [@(1).val.children()]),
     DirectSum_Single := Rule([DirectSum, @(1)], e -> @(1).val),
     Compose_Assoc := ARule(Compose, [@(1,Compose)], e -> [@(1).val.children()]),
     Compose_Single := Rule([Compose, @(1)], e -> @(1).val),
     J_J := ARule(Compose, [@(1,J), @(2,J)], e -> [ I(@(1).val.params[1]) ]),

     # ===============
     # Diags
     # ================
     Diag_Diag := ARule(Compose, [[Diag, @(1)], [Diag, @(2)]], e ->
         [ Diag(diagMul(@(1).val, @(2).val)) ]),

     Diag_L_Diag := ARule(Compose, [[Diag, @(1)], @(2, L), [Diag, @(3)]], e ->
         [ Diag(diagMul(@(1).val, fCompose(@(3).val, @(2).val))) * @(2).val ]),

     Tensor_Diag := Rule([Tensor, [I, @(1)], [Diag, @(2)]],
         e -> Diag(diagTensor(fConst(@(1).val, 1), @(2).val))),
     DirectSum_Diag := Rule([DirectSum, [I, @(1)], [Diag, @(2)]],
         e -> Diag(diagDirsum(fConst(@(1).val, 1), @(2).val))),

     E1_Diag_E1 := ARule(Compose, [[@(1,GathExtend,e->e.transposed), @, Odd1],
                                   [Diag, @(2)],
                                   [@(1,GathExtend,e->not e.transposed), @, Odd1]],
      e -> let(n:=@(2).val.domain(), f := @(2).val,
           [Diag(diagMul(fConst(n/2, 1/4),
                       diagAdd(fCompose(f, fTensor(fBase(2, 0), fId(n/2))),
                               fCompose(f, fTensor(fBase(2, 1), J(n/2))))))])),

     Scale_Scale := Rule([Scale, @(1), [Scale, @(2), @(3)]], e->Scale(@(1).val * @(2).val, @(3))),
     Scale_1 := Rule([Scale, @(1).cond(e->e=1), @(2)], e->@(2).val),

     # =============
     # PushL/PushR
     # =============
     PushL_within_PushL := Rule(@@(1,PushL, (e,cx)->IsBound(cx.PushL) and cx.PushL<>[]), e -> e.child(1)),

     Compose_PushL_PushL := ARule(Compose, [[PushL, @(1)], [PushL, @(2)]], e -> [ PushL(@(1).val * @(2).val) ]),
     Compose_PushR_PushL := ARule(Compose, [[PushR, @(1)], [PushL, @(2)]], e -> [ PushR(@(1).val * @(2).val) ]),
     XXX_PushL := ARule(Compose, [@(1, [Diag, L]), [PushL, @(2)]], e -> [PushL(@(1).val * @(2).val)]),

     Tensor_PushL := Rule([Tensor, @(1, I), [Compose, @(2, PushL), @(3)]],
	 e -> PushL(Tensor(@(1).val, @(2).val.child(1))) * Tensor(@(1).val, @(3).val)),

     DirectSum_PushL := Rule(@(1, DirectSum,
	     e -> ForAny(e.children(),
		 c -> ObjId(c)=PushL or (ObjId(c)=Compose and ObjId(c.child(1)) = PushL))),

	 e -> let(ch := @(1).val.children(),
	          split := List(ch, e -> Cond(ObjId(e) = PushL,
			                          [e, I(Cols(e))],
			                      ObjId(e) = Compose and ObjId(e.child(1))=PushL,
					          [e.child(1).child(1), Drop(e.children(), 1)],
					      [I(Rows(e)), e])),
		  PushL(DirectSum(List(split, x->x[1]))) * DirectSum(List(split, x->x[2]))))
));

RewriteRules(RulesDFTSymmetry, rec(
     BRDFT3_ExtendOdd0 := ARule(Compose, [[BRDFT3, @(1).cond(e->IsEvenInt(e)), 1/4], [GathExtend, @, Odd0]],
	 e -> let(k:=@(1).val,
	     [ GathExtend(k,Upsample0),SkewDTT(DCT3(k/2),1/2)])),

     Odd1_split := ARule(Compose, [[Tensor, [I,@(1)],@(2)],[GathExtend, @, Odd1]],
	 e -> let(n:=@(1).val,A:=@(2).val,k:=When(IsEvenInt(n),n/2, (n-1)/2),When(IsEvenInt(n),
	     [VStack(Tensor(I(k),A)*(1/2),Tensor(J(k),A*J(Cols(A))*(-1/2)))],
	     [VStack(DirectSum(Tensor(I(k),A)*(1/2),A*GathExtend(Cols(A),Odd1)),Tensor(J(k),A*J(Cols(A))*(-1/2)))]))),

     ReduceUpsample0_RC := ARule(Compose, [[ScatReduce, @(1), Upsample0],[RC, @(2)]],
	 e -> [@(2).val*ScatReduce(@(1).val,Upsample0)]),

     Push_ScatReduceWithinSum :=  ARule(Compose, [[ScatReduce, @, Upsample0],[@(1,IterDirectSum).cond(e->IsEvenInt(e.domain)),@(2)]],
	 e -> [IterDirectSum(@(1).val.var,@(1).val.domain,ScatReduce(Rows(@(2).val),Upsample0)*@(2).val)]
),

     Cut_SkewPRDFT := ARule(Compose, [[ScatReduce, @, Upsample0],[BSkewPRDFT,@(1),@(2)]],
	 e ->  let(n:=@(1).val,a:=@(2).val,f:=2*cospi(2*a),
	        [SkewDTT(DCT3(n/2),2*a)*When(IsEvenInt(n/2),
		     HStack(I(n/2),DirectSum(Mat([[f/2]]),Tensor(I(((n/2)-2)/2),Mat([[f,-1],[-1,f]])),Mat([[f-1]]))^M((n)/2,n/4)),
		     HStack(I(n/2),DirectSum(Mat([[f/2]]),Conjugate(Tensor(I((n/2-1)/2),Mat([[f,-1],[-1,f]])),M(n/2-1,(n/2-1)/2)))))])),

     SplitIterDirectSum := Rule([@(1,IterDirectSum), [Compose,@(2,SkewDTT),@(3,HStack)]],
	 e->let(i:=@(1).val.var,IterDirectSum(@(1).val.var,@(1).val.domain,@(2).val)*IterDirectSum(@(1).val.var,@(1).val.domain,@(3).val))),

     BRDFT3_ExtendOdd1 := ARule(Compose, [[BRDFT3, @(1).cond(e->IsEvenInt(e)), 1/4], [GathExtend, @, Odd1]],
         e -> let(k:=@(1).val,i:=Ind(),
	        [IterDirectSum(i,k/2,Mat([[1/2+cospi((2*i+1)/k)],[-1/2]]))*PolyDTT(DCT4(k/2))])),

     BRDFT3_J := ARule(Compose, [[BRDFT3, @(1).cond(e->IsEvenInt(e)), 1/4], [J, @]],
         e -> let(k:=@(1).val,i:=Ind(),
                [IterDirectSum(i,k/2,Mat([[-2*cospi((2*i+1)/k),-1-2*cospi((2*i+1)/(k/2))],[1,-2*cospi((2*((k/2)-i)-1)/k)]]))*BRDFT3(k)]))
));

#z := L(15,3) * GathExtend(15, Odd0);; zz := RulesDFTSymmetry(z);; PrintMat(MatSPL(z)-MatSPL(zz));
#z := L(16,4) * GathExtend(16, Odd0);; zz := RulesDFTSymmetry(z);; PrintMat(MatSPL(z)-MatSPL(zz));

#z := L(15,3) * GathExtend(15, Even00);; zz := RulesDFTSymmetry(z);; PrintMat(MatSPL(z)-MatSPL(zz));
#z := L(16,4) * GathExtend(16, Odd1);; zz := RulesDFTSymmetry(z);; PrintMat(MatSPL(z)-MatSPL(zz));

#RulesDFTSymmetry(Tensor(I(3), DFT(5))*L(15,3)*GathExtend(15, Even0));
#RulesDFTSymmetry(Tensor(I(2), DFT3(8))*L(16,2)*GathExtend(16, Even0));

dft1 := (N, k) -> Tensor(DFT1(k), I(N/k)) * Diag(Tw1(N, N/k, 1)) * Tensor(I(k), DFT1(N/k)) * L(N, k);
dft2 := (N, k) -> Tensor(DFT2(k), I(N/k)) * Diag(Tw2(N, N/k, 1)) * Tensor(I(k), DFT1(N/k)) * L(N, k);
dft3 := (N, k) -> Tensor(DFT1(k), I(N/k)) * Diag(Tw3(N, N/k, 1)) * Tensor(I(k), DFT3(N/k)) * L(N, k);
dft4 := (N, k) -> Tensor(DFT2(k), I(N/k)) * Diag(Tw4(N, N/k, 1)) * Tensor(I(k), DFT3(N/k)) * L(N, k);

pdft1 := (N, k, lsym, rsym) -> GathExtend(N, lsym).transpose() * dft1(N, k) * GathExtend(N, rsym);
pdft2 := (N, k, lsym, rsym) -> GathExtend(N, lsym).transpose() * dft2(N, k) * GathExtend(N, rsym);
pdft3 := (N, k, lsym, rsym) -> GathExtend(N, lsym).transpose() * dft3(N, k) * GathExtend(N, rsym);
pdft4 := (N, k, lsym, rsym) -> GathExtend(N, lsym).transpose() * dft4(N, k) * GathExtend(N, rsym);

# DCT6(8)
f := L(15, 3) * Tensor(I(5), DFT2(3, 1)) * L(15, 5) *
     Diag(Tw2(15, 5, 1)) *
     Tensor(I(3), DFT(5, 1)) * L(15, 3);
ff := GathExtend(Rows(f), Odd0).transpose() * f * GathExtend(Cols(f), Even1);

# f := L(16,4) * Tensor(I(4), DFT(4)) * L(16,4) * Diag(Tw1(16,4,1)) * Tensor(I(4), DFT(4)) * L(16,4);
# ff := GathExtend(Rows(f), Odd00).transpose() * f * GathExtend(Cols(f), Odd00);

# DCT2/DCT6
# f := L(16,4) * Tensor(I(4), DFT2(4)) * L(16,4) * Diag(Tw2(16,4)) * Tensor(I(4), DFT(4)) * L(16,4);
#ff := GathExtend(Rows(f), Odd0).transpose() * f * GathExtend(Cols(f), Even1);

# f4 := L(16,4) * Tensor(I(4), DFT2(4)) * L(16,4) * Diag(Tw4(16,4)) * Tensor(I(4), DFT3(4)) * L(16,4);
# ff := GathExtend(Rows(f), Odd1).transpose() * f4 * GathExtend(Cols(f), Odd1);

#  fo := L(21,3) * Tensor(I(7), DFT2(3)) * L(21,7) * Diag(Tw4(21,7)) * Tensor(I(3), DFT3(7)) * L(21,3);
# ffo := GathExtend(Rows(fo), Odd1).transpose() * fo * GathExtend(Cols(fo), Odd1);

# fr1 := GathExtend(16, CE0_R2Cpx(W_RFT, 1).invertW()).transpose() * f * GathExtend(16, Real);

 fr2 := GathExtend(21, CE0_R2Cpx(W_RFT, 1).invertW()).transpose() *
        L(21,3) * Tensor(I(7), DFT(3)) * L(21,7) * Diag(Tw1(21,7,1)) * Tensor(I(3), DFT(7)) * L(21,3) *
	GathExtend(21, Real);

 fr4 := GathExtend(21, CO1_R2Cpx(W_RFT, -E(4)).invertW()).transpose() *
        L(21,3) * Tensor(I(7), DFT2(3)) * L(21,7) * Diag(Tw4(21,7,1)) * Tensor(I(3), DFT3(7)) * L(21,3) *
	GathExtend(21, Real);

# fr4 := GathExtend(16, CO1_R2Cpx).transpose() * L(16,4) * Tensor(I(4), DFT2(4)) * L(16,4) * Diag(Tw4(16,4)) * Tensor(I(4), DFT3(4)) * L(16,4) * GathExtend(16, Real);

#p  := Mat(1/2*[[1,  E(4)]]);
#pp := Mat(1/2*[[1, -E(4)]]);

#them := MatSPL(DFT(4) * Tensor(I(4), 2*pp));;
#me   := MatSPL(Tensor(I(4), 2*pp) * RC(DFT(4,1)));;

# y = DFT x* = (DFT* x)* ?
#
# C IxP = IxP' RC(C)
# IxP' = C IxP RC(C)^-1

# C IxP =IxP' C Ix(1,j)

dftproj := function(f)
    f := RulesDFTSymmetry(f);
    f := f.transpose();
    f := RulesDFTSymmetry(f);
    f := f.transpose();
    return f;
end;


text:
#############################################################################
####
##
#A  anustab.gi                  ANUPQ package                  Eamonn O'Brien
#A                                                              Werner Nickel
##
#Y  Copyright 1993-2001,  Lehrstuhl D fuer Mathematik,  RWTH Aachen,  Germany
#Y  Copyright 1993-2001,  School of Mathematical Sciences, ANU,     Australia
##
#W  Greg Gamble reformulated the original code as a function and  then  split
#W  the original `anustab.g' into the declare/install files anustab.g[di].
##
##  Install file for function to  compute  the  stabiliser  of  an  allowable
##  subgroup; description is written to file LINK_output.
##

#############################################################################
##
#F  PqStabiliserOfAllowableSubgroup( <ANUPQglb>, <F>, <gens>, <relativeOrders>,
##  <ANUPQsize>, <ANUPQagsize>) . .  compute stabiliser of allowable subgroup
##
##  computes stabiliser and writes `LINK_output' output  file  for  the  `pq'
##  binary. The `LINK_output' file is written to `ANUPQData.tmpdir', but when
##  the `pq' binary was started up directly by the user that directory is set
##  to    the    current    directory.    The    argument    variables     of
##  `PqStabiliserOfAllowableSubgroup' are essentially the variables  assigned
##  to by `GAP_link_via_file.c' in the `pq' source code.
##
InstallGlobalFunction( PqStabiliserOfAllowableSubgroup,
function( ANUPQglb, F, gens, relativeOrders, ANUPQsize, ANUPQagsize )
local m, n, H, pcgs, aut, NumberAgAutos, i, imgs, NumberGlAutos, p, d,
      ANUPQMaxDegree, V, elm, baseU, id, baseN, infoLevelAutGrp, LINK_output,
      soluble, a, mat;

  # number of generators
  if Length( ANUPQglb.glAutos ) > 0 then
      m := Length (ANUPQglb.glAutos[1]);
  elif Length( ANUPQglb.agAutos ) > 0 then
      m := Length (ANUPQglb.agAutos[1]);
  else
      Error( "Panic: PqStabiliserOfAllowableSubgroup ",
             "called without automorphisms" );
  fi;

  # get the p-group
  n := Length (Pcgs(F));
  H := F / Subgroup (F, Pcgs(F){[m+1 .. n]}); IsPGroup(H);
  pcgs := Pcgs(H);

  # set up automorphism rec
  aut := rec();

  NumberAgAutos := Length (ANUPQglb.agAutos);
  aut.agAutos := [1..NumberAgAutos];
  for i in [1..Length (ANUPQglb.agAutos)] do
    MakeImmutable( ANUPQglb.genQ[i] );
    ConvertToMatrixRep( ANUPQglb.genQ[i], ANUPQglb.F );
    imgs := List( ANUPQglb.agAutos[i], x -> PcElementByExponents( pcgs, x ) );
    aut.agAutos[i] := PGAutomorphism( H, pcgs, imgs);
    aut.agAutos[i]!.mat := ANUPQglb.genQ[i];
  od;
  aut.agAutos := Reversed( aut.agAutos );
  aut.agOrder := Reversed (relativeOrders);

  NumberGlAutos := Length( ANUPQglb.glAutos );
  aut.glAutos := [1..NumberGlAutos];
  for i in [1..NumberGlAutos] do
    MakeImmutable( ANUPQglb.genQ[i] );
    ConvertToMatrixRep( ANUPQglb.genQ[NumberAgAutos + i], ANUPQglb.F );
    imgs := List( ANUPQglb.glAutos[i], x -> PcElementByExponents( pcgs, x ) );
    aut.glAutos[i] := PGAutomorphism( H, pcgs, imgs);
    aut.glAutos[i]!.mat := ANUPQglb.genQ[NumberAgAutos+ i];
  od;

  aut.field    := ANUPQglb.F;
  aut.group    := H;
  if ANUPQsize <> fail then
    aut.size     := ANUPQsize;
  fi;
  aut.one      := IdentityPGAutomorphism( H );
  aut.one!.mat := IdentityMat( ANUPQglb.q, aut.field );

  # add perm oper
  p := Size( aut.field );
  d := RankPGroup( aut.group );
  #compute perm rep having at most this degree
  ANUPQMaxDegree := 10000;
  if p^d <= ANUPQMaxDegree or ANUPQsize = fail then
    V := aut.field^d;
    if ANUPQsize = fail and p^d > ANUPQMaxDegree then
        Info(InfoANUPQ + InfoWarning, 1,
             "no. of elements in vector space V is large (", p^d, ")");
        Info(InfoANUPQ + InfoWarning, 1,
             "... this may take a while, if it succeeds at all!");
    fi;
    elm := AsSet( V );
    aut.glOper := [];
    for i in [1..NumberGlAutos] do
        a := aut.glAutos[i];
        mat := List(a!.baseimgs,
                    x -> ExponentsOfPcElement( pcgs, x ){[1..d]});
        mat := mat * One( aut.field );
        MakeImmutable( mat );
        ConvertToMatrixRep( mat, aut.field );
        aut.glOper[i] := Permutation( mat, elm, OnRight );
    od;
    #PrintTo("perms", "aut.glOper := ", aut.glOper,"; \n" );
    aut.glOrder := Size( Group( aut.glOper, () ) );
  else
    aut.glOrder := ANUPQsize / Size(ANUPQglb.F)^(ANUPQagsize);
  fi;

  Info(InfoANUPQ, 2, "Order of GL subgroup is ", aut.glOrder);
  Info(InfoANUPQ, 2, "No. of soluble autos is ", Length (aut.agAutos));

  # get allowable subgroup and nucleus
  baseU := gens * One (aut.field);
  id := IdentityMat( ANUPQglb.q, aut.field);
  baseN := id{[1..ANUPQglb.r]};

  # call stabilizer
  infoLevelAutGrp := InfoLevel( InfoAutGrp );
  if InfoLevel( InfoANUPQ ) >= 2 then
    SetInfoLevel( InfoAutGrp, 3 );
  fi;
  PGOrbitStabilizer( aut, baseU, baseN, false );

  if (Length (aut.glAutos) = 0) then
    soluble := -1;
  else
    soluble := 0;
  fi;

  aut.agAutos := Reversed (aut.agAutos);
  aut.agOrder := Reversed (aut.agOrder);

  LINK_output := Filename( ANUPQData.tmpdir, "LINK_output" );
  PrintTo (LINK_output, soluble);
  AppendTo (LINK_output, "\n");
  AppendTo (LINK_output, Length (aut.agAutos));
  AppendTo (LINK_output, "\n");
  for i in [1..Length (aut.agOrder)] do
    AppendTo (LINK_output, aut.agOrder[i]);
    AppendTo (LINK_output, "\n");
  od;
  AppendTo (LINK_output, Length (aut.agAutos) + Length (aut.glAutos) );
  AppendTo (LINK_output, "\n");
  for a in aut.agAutos do
    mat := List( a!.baseimgs, x -> ExponentsOfPcElement( pcgs, x ) );
    mat := Flat( mat );
    for m in mat do
        AppendTo (LINK_output, m, " " );
    od;
    AppendTo (LINK_output, "\n");
  od;
  for a in aut.glAutos do
    mat := List( a!.baseimgs, x -> ExponentsOfPcElement( pcgs, x ) );
    mat := Flat( mat );
    for m in mat do
        AppendTo (LINK_output, m, " " );
    od;
    AppendTo (LINK_output, "\n");
  od;
  SetInfoLevel( InfoAutGrp, infoLevelAutGrp );
end );

#E  anustab.gi  . . . . . . . . . . . . . . . . . . . . . . . . . . ends here

ext: gi

max_stars_repo_name: gap-system/anupq

max_issues_repo_name: gap-system/anupq

max_forks_repo_name: gap-system/anupq


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

# NOTE: find a better place for these! What about a single standard drop tag rule?
NewRulesFor(TRC, rec(
    TRC_tag := rec(
        forTransposition := false,
        applicable := (self, nt) >>
            (nt.isTag(1, spiral.paradigms.smp.AParSMP) or not nt.hasTags())
            # AVecReg is taken from a namespace that is NOT YET LOADED
	    # hence the fully qualified name
            and not nt.hasTag(spiral.paradigms.vector.AVecReg)
            and not nt.hasTag(spiral.paradigms.vector.AVecRegCx),

        children := nt -> [[ nt.params[1].withTags(nt.getTags()) ]],
        apply := (nt, c, cnt) -> RC(c[1])
    )
));

NewRulesFor(TDiag, rec(
    TDiag_tag := rec(
        forTransposition := false,

	# YSV: Below limits applicability to the cases where diag size is divisible by vlen
	#      which is a safe thing to do. Because VectorCodegen, can't  generated code
	#      for VDiags of size non-divisible by vlen. HOWEVER, if VDiag is propagate
	#      past any kind of VGath, this problem goes away. So having no restriction,
	#      will work MOST of the time, but not all the time.
        #
	# applicable := (self, nt) >> let(
	#     vtags := [spiral.paradigms.vector.AVecReg, spiral.paradigms.vector.AVecRegCx],
	#     dom   := nt.params[1].domain(),
	#     not nt.hasAnyTag(vtags) or (dom mod nt.getAnyTag(vtags).v) = 0
	# ),

        apply := (t, C, Nonterms) -> let(
	    vtags := [spiral.paradigms.vector.AVecReg, spiral.paradigms.vector.AVecRegCx],
	    Cond(t.hasAnyTag(vtags),
		 spiral.paradigms.vector.sigmaspl.VDiag(t.params[1], t.getAnyTag(vtags).v),
		 Diag(t.params[1])
	    )
	)
    )
));

RulesFor(TRCDiag, rec(
    TRCDiag_tag := rec(
        forTransposition := false,
        applicable := (self, nt) >> not nt.transposed,
        rule := (P, C) -> RC(Diag(P[1])))
));

RulesFor(TId, rec(
    TId_tag := rec(
        forTransposition := false,
        switch := false,
        rule := (P, C) -> P[1])
));

NewRulesFor(TRaderMid, rec(
    TRaderMid_tag := rec(
        forTransposition := false,
        apply := (t, C, Nonterms) -> t.raderMid(t.params[1], t.params[2], t.params[3])
    )
));

NewRulesFor(TRDiag, rec(
    TRDiag_RT_Diag := rec(
        forTransposition := true,
        apply := (t, C, Nonterms) -> t.terminate()
    )
));

NewRulesFor(TCompose, rec(
    TCompose_tag := rec(
        forTransposition := false,
        applicable := (self, nt) >> true,
        children := nt -> [ List(nt.params[1], e -> e.withTags(nt.getTags())) ],
        apply := (nt, c, nt) -> Grp(Compose(c))
    )
));

NewRulesFor(TCond, rec(
    TCond_tag := rec(
        forTransposition := false,
        applicable := (self, nt) >> true,
        children := nt -> [[
	    nt.params[2].withTags(nt.getTags()), nt.params[3].withTags(nt.getTags()) ]],
        apply := (t, C, Nonterms) -> COND(t.params[1], C[1], C[2])
    )
));

NewRulesFor(TGrp, rec(
    TGrp_tag := rec(
        forTransposition := false,
        children := nt -> [[ nt.params[1].withTags(nt.getTags()) ]],
        apply := (nt, c, cnt) -> Grp(c[1])
    )
));

NewRulesFor(TInplace, rec(
    TInplace_tag := rec(
        forTransposition := false,
        children := nt -> [[ nt.params[1].withTags(nt.getTags()) ]],
        apply := (nt, c, cnt) -> Inplace(c[1])
    )
));

NewRulesFor(TICompose, rec(
    TICompose_tag := rec(
        forTransposition := false,
        children := nt -> [[ nt.params[3].withTags(nt.getTags()) ]],
        apply := (nt, c, cnt) -> ICompose(nt.params[1], nt.params[2], c[1])
    )
));

########################################################################
#   (A + B) rules
NewRulesFor(TDirectSum, rec(
#   (A + B) terminate
    A_dirsum_B := rec(
        forTransposition := false,
        children := (self, t) >> let( tags := t.getTags(),
            [[ t.params[1].withTags(tags), t.params[2].setTags(tags) ]]
        ),
        apply := (t, C, Nonterms) -> DirectSum(C)

#D        children := (self, t) >> let (tags:=GetTags(t),
#D            [[ AddTag(t.params[1], tags), SetTag(t.params[2], tags) ]]),
    )
));

########################################################################
#   (A x B) rules
NewRulesFor(TTensor, rec(
#   (A x B) -> (A x I)(I x B)
    AxI_IxB := rec(
        info := "(A x B) -> (A x I)(I x B)",
        forTransposition := false,
        applicable := nt -> true,
        inplace := false,
        children := (self, nt) >> let(inp := When(self.inplace, TInplace, x->x),
            [[ TCompose([
                inp(TTensorI(nt.params[1], nt.params[2].dims()[1], AVec, AVec)),
                TTensorI(nt.params[2], nt.params[1].dims()[2], APar, APar)
            ]).withTags(nt.getTags()) ]]),
        apply := (nt, c, cnt) -> c[1],
#D        isApplicable := P -> true,
#D        allChildren := P -> [[TCompose([TTensorI(P[1], P[2].dims()[1], AVec, AVec), TTensorI(P[2], P[1].dims()[2], APar, APar)], P[3])]],
#D        rule := (P, C) -> C[1]
    ),
#   (A x B) -> (I x B)(A x I)
    IxB_AxI := rec(
        info := "(A x B) -> (I x B)(A x I)",
        forTransposition := false,
        applicable := nt -> true,
        inplace := false,
        children := (self, nt) >> let(inp := When(self.inplace, TInplace, x->x),
            [[ TCompose([
                inp(TTensorI(nt.params[2], nt.params[1].dims()[1], APar, APar)),
                TTensorI(nt.params[1], nt.params[2].dims()[2], AVec, AVec)
            ]).withTags(nt.getTags()) ]]),
        apply := (nt, c, cnt) -> c[1]

#D        isApplicable := P -> true,
#D        allChildren := P -> [[TCompose([TTensorI(P[2], P[1].dims()[1], APar, APar), TTensorI(P[1], P[2].dims()[2], AVec, AVec)], P[3])]],
#D        rule := (P, C) -> C[1]
    ),
#   (A x B) -> (L(B x I))(L(A x I))
    L_BxI__L_AxI := rec(
        info := "(A x B) -> (L(B x I))(L(A x I))",
        forTransposition := false,
        applicable := nt -> true,
        children := nt -> [[ TCompose([
            TTensorI(nt.params[2], nt.params[1].dims()[1], APar, AVec),
            TTensorI(nt.params[1], nt.params[2].dims()[2], APar, AVec)
        ]).withTags(nt.getTags()) ]],
        apply := (nt, c, cnt) -> c[1]

#D        isApplicable := P -> true,
#D        allChildren := P -> [[TCompose([TTensorI(P[2], P[1].dims()[1], APar, AVec), TTensorI(P[1], P[2].dims()[2], APar, AVec)], P[3])]],
#D        rule := (P, C) -> C[1]
    ),
#   (A x B) -> ((A x I)L)((B x I)L)
    AxI_L__BxI_L := rec(
        info := "(A x B) -> ((A x I)L)((B x I)L)",
        forTransposition := false,
        applicable := nt -> true,
        children := nt -> [[ TCompose([
            TTensorI(nt.params[1], nt.params[2].dims()[1], AVec, APar),
            TTensorI(nt.params[2], nt.params[1].dims()[2], AVec, APar)
        ]).withTags(nt.getTags()) ]],
        apply := (nt, c, cnt) -> c[1],

#D        isApplicable := P -> true,
#D        allChildren := P -> [[TCompose([TTensorI(P[1], P[2].dims()[1], AVec, APar), TTensorI(P[2], P[1].dims()[2], AVec, APar)], P[3])]],
#D        rule := (P, C) -> C[1]
    ),
));

########################################################################
#   rules for A x I, I x A, (A x I)L, (I x A)L
NewRulesFor(TTensorI, rec(
    TTensorI_toGT := rec(
        applicable := t -> true,
        freedoms := t -> [], # no degrees of freedom
        child := (t, fr) -> [ GT_TTensorI(t) ], # fr will be an empty list
        apply := (t, C, Nonterms) -> C[1]
    )
));

NewRulesFor(TTensorI, rec(
#   base cases
#   I x A
    IxA_base := rec(
        info := "IxA base",
        forTransposition := false,
        applicable := nt -> (not nt.hasTags() or nt.firstTag() = ANoTag) and IsParPar(nt.params),
        children := nt -> [[ nt.params[1].withTags(nt.getTags()) ]],
        apply := (nt, c, cnt) -> When(nt.params[2] > 1,
            Tensor(I(nt.params[2]), c[1]),
            c[1]
        )
#D        isApplicable := (self, P) >> PUntagged(self.nonTerminal, P) and IsParPar(P),
#D        allChildren := P -> [[P[1]]],
#D        rule := (P, C) -> When(P[2]>1,Tensor(I(P[2]),C[1]),C[1])
    ),
#   A x I
    AxI_base := rec(
        info := "AxI base",
        forTransposition := false,
        applicable := nt -> (not nt.hasTags() or nt.firstTag() = ANoTag) and IsVecVec(nt.params),
        children := nt -> [[ nt.params[1].withTags(nt.getTags()) ]],
        apply := (nt, c, cnt) -> When( nt.params[2] > 1,
            Tensor(c[1], I(nt.params[2])),
            c[1]
        ),
#D        isApplicable := (self, P) >> PUntagged(self.nonTerminal, P) and IsVecVec(P),
#D        allChildren := P -> [[P[1]]],
#D        rule := (P, C) -> When(P[2]>1,Tensor(C[1], I(P[2])),C[1])
    ),
#   (I x A)L
    IxA_L_base := rec(
        info := "(IxA)L base",
        forTransposition := false,
        applicable := nt -> (not nt.hasTags() or nt.firstTag() = ANoTag) and IsParVec(nt.params),
        children := nt -> [[ nt.params[1].withTags(nt.getTags()) ]],
        apply := (nt, c, cnt) -> Tensor(I(nt.params[2]), c[1]) * L(c[1].dims()[2] * nt.params[2], nt.params[2]),

#D        isApplicable := (self, P) >> PUntagged(self.nonTerminal, P) and IsParVec(P),
#D        allChildren := P -> [[P[1]]],
#D        rule := (P, C) -> Tensor(I(P[2]), C[1])*L(C[1].dims()[2]*P[2], P[2])
    ),
#   L(I x A)
    L_IxA_base := rec(
        info := "L(IxA) base",
        forTransposition := false,
        applicable := nt -> (not nt.hasTags() or nt.firstTag() = ANoTag) and IsVecPar(nt.params),
        children := nt -> [[ nt.params[1].withTags(nt.getTags()) ]],
        apply := (nt, c, cnt) -> L(c[1].dims()[1] * nt.params[2], c[1].dims()[1]) * Tensor(I(nt.params[2]), c[1])

#D        isApplicable := (self, P) >> PUntagged(self.nonTerminal, P) and IsVecPar(P),
#D        allChildren := P -> [[P[1]]],
#D        rule := (P, C) -> L(C[1].dims()[1]*P[2], C[1].dims()[1]) * Tensor(I(P[2]), C[1])
    ),
#   splitting rules ##############################################################
#   (A _m x I_n)L_mn_m
    AxI_L_split := rec(
        info := "split (A_m x I_n) L^mn_m --> (L_mn/u_m x I_u) * (I_n/u x (A_m x I_u) * L_mu_m )",
        forTransposition := false,
        applicable := nt -> (nt.firstTag().kind() = AGenericTag) and IsVecPar(nt.params),
        children := nt -> let(t := nt.getTags(), p := nt.params, d  := p[1].dims(), mu := t[1].params[1], [
            TTensorI(TL(d[1] * p[2]/mu, d[1],1,1), mu, AVec, AVec).withTags(t),
            TTensorI(p[1], mu, AVec, APar).withTags(t)
        ]),
        apply := (nt, c, cnt) -> let(t := nt.getTags(), n := nt.params[2], mu := t[1].params[1],
            c[1] * Tensor(I(n/mu), c[2])
        ),

        # Example
        # =======
        # t:=TTensorI(DFT(4, 1), 4, AVec, APar).withTags([ AGenericTag(2) ]);
        # c:=AxI_L_split.children(t);
        # res := AxI_L_split.apply(t,c,false);

        switch:=false
    ),

#   (I_n x A_rxs) L^ns_n
    IxA_L_split := rec(
        info := "split (I_n x A_rxs) L^ns_n",
        forTransposition := false,
        applicable := nt -> IsParVec(nt.params),
        children := nt -> let(t := nt.getTags(), p := nt.params, d := p[1].dims(), [[
            TTensorI(p[1], p[2], APar, APar).withTags(t),
            TL(d[2]*p[2], p[2], 1, 1).withTags(t)
        ]]),
        apply := (nt, c, cnt) -> c[1] * c[2],

#D        isApplicable := P -> P[3].isPar and P[4].isVec,
#D        allChildren := P -> let(pv:=P[5], d:=P[1].dims(), [[TTensorI(P[1], P[2], APar, APar, pv), TL(d[2]*P[2], P[2], 1, 1, pv)]]),
#D        rule := (P, C) -> C[1] * C[2],
        switch := false
    ),
#   L^nr_n (A_rxs x I_n)
    L_AxI_split := rec(
        info := "split L^nr_n (A_rxs x I_n) ",
        forTransposition := false,
        applicable := nt -> IsParVec(nt.params),
        children := nt -> let( t := nt.getTags(), p := nt.params, d := p[1].dims(), [[
            TL(d[1] * p[2], p[2], 1, 1).withTags(t),
            TTensorI(p[1], p[2], AVec, AVec).withTags(t)
        ]]),
        apply := (nt, c, cnt) -> c[1] * c[2],
        switch := false
#D        isApplicable := P -> P[3].isPar and P[4].isVec,
#D        allChildren := P -> let(pv:=P[5], d:=P[1].dims(), [[ TL(d[1]*P[2], P[2], 1, 1, pv), TTensorI(P[1], P[2], AVec, AVec, pv) ]]),
#D        rule := (P, C) -> C[1] * C[2],
    ),
#   L^nr_r (I_n x A_rxs)
    L_IxA_split := rec(
        info := "split L^nr_r (I_n x A_rxs)",
        forTransposition := false,
        applicable := nt -> IsVecPar(nt.params),
        children := nt -> let( t := nt.getTags(), p := nt.params, d := p[1].dims(), [[
            TL(d[1]*p[2], d[1], 1, 1).withTags(t),
            TTensorI(p[1], p[2], APar, APar).withTags(t)
        ]]),
        apply := (nt, c, cnt) -> c[1] * c[2],

#D        isApplicable := P -> P[3].isVec and P[4].isPar,
#D        allChildren := P -> let(pv:=P[5], d:=P[1].dims(), [[TL(d[1]*P[2], d[1], 1, 1, pv), TTensorI(P[1], P[2], APar, APar, pv)]]),
#D        rule := (P, C) -> C[1] * C[2],
        switch := false
    ),
#   (A_rxs x I_n) L^nr_s
    AxI_L_split := rec(
        info := "split (A_rxs x I_n) L^nr_s ",
        forTransposition := false,
        applicable := nt -> IsVecPar(nt.params),
        children := nt -> let( t := nt.getTags(), p := nt.params, d := p[1].dims(), [[
            TTensorI(p[1], p[2], APar, APar).withTags(t),
            TL(d[2]*p[2], d[2], 1, 1).withTags(t)
        ]]),
        apply := (nt, c, cnt) -> c[1] * c[2],
#D        isApplicable := P -> P[3].isVec and P[4].isPar,
#D        allChildren := P -> let(pv:=P[5], d:=P[1].dims(), [[ TTensorI(P[1], P[2], APar, APar, pv), TL(d[2]*P[2], d[2], 1, 1, pv)]]),
#D        rule := (P, C) -> C[1] * C[2],
        switch := false
    ),
##   vector recursion #############################################################
#   (I x (I x A)L)L
    IxA_L_vecrec := rec(
        info := "(I x (I x A)L)L vector recursion",
        forTransposition := false,
        applicable := nt -> ObjId(nt.params[1]) = TTensorI and IsParVec(nt.params) and IsParVec(nt.params[1].params),
        children := nt -> let(k := nt.params[2], m := nt.params[1].params[2], n := nt.params[1].params[1].dims(), [[
            TL(k*m, k, 1, n[1]).withTags(nt.getTags()),
            TTensorI(nt.params[1].params[1], nt.params[2], APar, AVec).withTags(nt.getTags()),
            TL(m*n[2], m, 1, k).withTags(nt.getTags())
        ]]),
        apply := (nt, c, cnt) -> let(m := nt.params[1].params[2],
            c[1] * Tensor(I(m), c[2]) * c[3]
        ),
#D        isApplicable := P -> P[1].name = "TTensorI" and P[3].isPar and P[4].isVec and P[1].params[3].isPar and P[1].params[4].isVec,
#D        allChildren := P -> let(k:=P[2], m:=P[1].params[2], n:=P[1].params[1].dims(),
#D                [[ TL(k*m, k, 1, n[1], P[5]), TTensorI(P[1].params[1], P[2], APar, AVec, P[5]),  TL(m*n[2], m, 1, k, P[5])]]),
#D        rule := (P, C) -> let(k:=P[2], m:=P[1].params[2], n:=P[1].params[1].dims(),
#D                C[1] * Tensor(I(m), C[2]) * C[3]
#D            ),
        switch := false
    ),
#   L(I x L(I x A))
    L_IxA_vecrec := rec(
        info := "L(I x L(I x A)) vector recursion",
        forTransposition := false,
        applicable := nt -> ObjId(nt.params[1]) = TTensorI and IsVecPar(nt.params) and IsVecPar(nt.params[1].params),
        children := nt -> let( k := nt.params[2], m := nt.params[1].params[2], n := nt.params[1].params[1].dims(), [[
            TL(m*n[1], n[1], 1, k).withTags(nt.getTags()),
            TTensorI(nt.params[1].params[1], nt.params[2], AVec, APar).withTags(nt.getTags()),
            TL(k*m, m, 1, n[2]).withTags(nt.getTags())
        ]]),
        apply := (nt, c, cnt) -> let(m := nt.params[1].params[2],
            c[1] * Tensor(I(m), c[2]) * c[3]
        ),
#D        isApplicable := P -> P[1].name = "TTensorI" and P[3].isVec and P[4].isPar and P[1].params[3].isVec and P[1].params[4].isPar,
#D        allChildren := P -> let(k:=P[2], m:=P[1].params[2], n:=P[1].params[1].dims(),
#D                [[ TL(m*n[1], n[1], 1, k, P[5]), TTensorI(P[1].params[1], P[2], AVec, APar, P[5]),  TL(k*m, m, 1, n[2], P[5])]]),
#D        rule := (P, C) -> let(k:=P[2], m:=P[1].params[2], n:=P[1].params[1].dims(),
#D                C[1] * Tensor(I(m), C[2]) * C[3]
#D            ),
        switch := false
    )
));

########################################################################
#   rules for L

#D isVec := P->Length(P[5]) > 0 and P[5][1].isVec;

NewRulesFor(TL, rec(
#   TL(N,n,l,r,[]) -> I_l x L(N,n) x I_r
    L_base := rec(
        forTransposition := false,
        applicable := nt -> nt.isTag(1, spiral.paradigms.smp.AParSMP) or not nt.hasTags(),
        apply := (nt, c, cnt) -> let(
            c1 := When(nt.params[3]=1, [], [I(nt.params[3])]),
            c2 := When(nt.params[4]=1, [], [I(nt.params[4])]),
            Tensor(Concat(c1, [ L(nt.params[1], nt.params[2]) ], c2))
        )
    ),
#   TL(N,n,l,r,[]) -> I_l x L(N,n) x I_r
    L_func := rec(
        forTransposition := false,
        applicable := nt -> nt.isTag(1, spiral.paradigms.smp.AParSMP) or not nt.hasTags(),
        apply := (nt, c, cnt) -> let(
            c1 := When(nt.params[3]=1, [], [fId(nt.params[3])]),
            c2 := When(nt.params[4]=1, [], [fId(nt.params[4])]),
            Prm(fTensor(Concat(c1, [ L(nt.params[1], nt.params[2]) ], c2)))
        )
    ),
#   recursion rules
    IxLxI_kmn_n := rec (
        info             := "I(l) x L(kmn, n) x I(r) -> (I_l x L(kn,n) x I(mr))(I(kl) x L(mn, n) x I(r))",
        forTransposition := false,
        applicable := nt -> Length(DivisorsIntDrop(nt.params[1]/nt.params[2])) > 0,
        children := nt -> let(
            N := nt.params[1], n := nt.params[2],
            km := N/n, ml := DivisorsIntDrop(km),
            l := nt.params[3], r := nt.params[4],
            List(ml, m -> let( k := km/m, [
                TL(k*n, n, l, r*m).withTags(nt.getTags()),
                TL(m*n, n, k*l, r).withTags(nt.getTags())
            ]))
        ),
        apply := (nt, c, cnt) -> let(
            spl := c[1] * c[2],
            When(nt.params[1] = nt.params[2]^2,
                SymSPL(spl),
                spl
            )
        ),

#D        isApplicable     := P -> #isVec(P) and let(v:=P[5][1].v, (P[1]*P[2] >= v or P[1]*P[3] >= v) and
#D                                Length(DivisorsIntDrop(P[1]/P[2])) > 0,
#D        allChildren := P -> let(N:=P[1], n:=P[2], km:=N/n, ml:=DivisorsIntDrop(km), l:=P[3], r:=P[4], vp:=P[5],
#D            List(ml, m->let(k:=km/m, [TL(k*n, n, l, r*m, vp), TL(m*n,n, k*l, r, vp)])) ),
#D        rule := (P, C) -> let(spl := C[1]*C[2], When(P[1]=P[2]^2, SymSPL(spl), spl)),
        switch := false
    ),
    IxLxI_kmn_km := rec (
        info             := "I(l) x L(kmn, km) x I(r) -> (I(kl) x L(mn,m) x I(r))(I(l) x L(kn, k) x I(r))",
        forTransposition := false,
        applicable := nt -> Length(DivisorsIntDrop(nt.params[2])) > 0,
        children := nt -> let(
            N := nt.params[1], km := nt.params[2],
            n := N/km, ml := DivisorsIntDrop(km),
            l := nt.params[3], r := nt.params[4],
            List(ml, m->let(
                k := km/m,
                [
                    TL(m*n, m, k*l, r).withTags(nt.getTags()),
                    TL(k*n,k, l, m*r).withTags(nt.getTags())
                ]
            ))
        ),
        apply := (nt, C, cnt) -> let(P := nt.params, spl := C[1]*C[2], When(P[1]=P[2]^2, SymSPL(spl), spl)),
#D        isApplicable     := P -> #isVec(P) and let(v:=P[5][1].v, (P[1]*P[2] >= v or P[1]*P[3] >= v) and
#D                                Length(DivisorsIntDrop(P[2])) > 0,
#D        allChildren := P -> let(N:=P[1], km:=P[2], n:=N/km, ml:=DivisorsIntDrop(km), l:=P[3], r:=P[4], vp:=P[5],
#D            List(ml, m->let(k:=km/m, [TL(m*n, m, k*l, r, vp), TL(k*n,k, l, m*r, vp)])) ),
#D        rule := (P, C) -> let(spl := C[1]*C[2], When(P[1]=P[2]^2, SymSPL(spl), spl)),
        switch := false
    ),
    IxLxI_IxLxI_up := rec (
        info             := "I(l) x L(kmn, km) x I(r) -> (I(l) x L(kmn, k) x I(r))(I(l) x L(kmn, m) x I(r))",
        forTransposition := false,
        applicable       := nt -> Length(DivisorPairs(nt.params[2])) > 0,
        children := nt -> let(
            N := nt.params[1], km := DivisorPairs(nt.params[2]),
            l := nt.params[3], r := nt.params[4], t := nt.getTags(),
            List(km, i->[TL(N, i[1], l, r).withTags(t), TL(N, i[2], l, r).withTags(t)])
        ),
        apply := (nt, c, nt) -> c[1] * c[2],

#D        isApplicable     := P -> Length(DivisorPairs(P[2])) > 0,
#D        allChildren := P -> let(N:=P[1], km:=DivisorPairs(P[2]), l:=P[3], r:=P[4], vp:=P[5],
#D            List(km, i->[TL(N, i[1], l, r, vp), TL(N, i[2], l, r, vp)])),
#D        rule := (P, C) -> C[1]*C[2],
        switch := false
    ),
    IxLxI_IxLxI_down := rec (
        info             := "I(l) x L(kmn, k) x I(r) -> (I(l) x L(kmn, km) x I(r))(I(l) x L(kmn, kn) x I(r))",
        forTransposition := false,
        applicable       := nt -> Length(DivisorPairs(nt.params[1]/nt.params[2])) > 0,
        children         := nt -> let(
            N := nt.params[1], km := DivisorPairs(nt.params[1]/nt.params[2]),
            l := nt.params[3], r := nt.params[4], t := nt.getTags(),
            List(km, i->[TL(N, N/i[1], l, r).withTags(t), TL(N, N/i[2], l, r).withTags(t)])
        ),
        apply := (nt, c, cnt) -> c[1] * c[2],

#D        isApplicable     := P -> Length(DivisorPairs(P[1]/P[2])) > 0,
#D        allChildren := P -> let(N:=P[1], km:=DivisorPairs(P[1]/P[2]), l:=P[3], r:=P[4], vp:=P[5],
#D            List(km, i->[TL(N, N/i[1], l, r, vp), TL(N, N/i[2], l, r, vp)])),
#D        rule := (P, C) -> C[1]*C[2],
        switch := false
    ),
    IxLxI_loop1 := rec(
        info := "I x L x I loop1",
        forTransposition := false,
        applicable := nt -> not nt.hasTags(),
        apply := (nt, c, cnt) -> let(
            m := nt.params[2], n := nt.params[1]/nt.params[2], j:=Ind(m), fid := fId(n), fbase := fBase(m,j),
            gath := Gath(fTensor(fid, fbase)), scat := Scat(fTensor(fbase, fid)),
            c0 := [ISum(j, m, scat*gath)],
            c1 := When(nt.params[3]=1, [], [I(nt.params[3])]),
            c2 := When(nt.params[4]=1, [], [I(nt.params[4])]),
            Tensor(Concat(c1,c0,c2))
        ),

#D        isApplicable := P -> Length(P[5]) = 0,
#D        rule := (P, C) -> let(m:=P[2], n:=P[1]/P[2], j:=Ind(m), fid := fId(n), fbase := fBase(m,j),
#D                gath := Gath(fTensor(fid, fbase)), scat := Scat(fTensor(fbase, fid)),
#D                C0 := [ISum(j, m, scat*gath)], C1:=When(P[3]=1, [], [I(P[3])]), C2:=When(P[4]=1, [], [I(P[4])]), Tensor(Concat(C1, C0, C2))),
        switch := false
    ),
    IxLxI_loop2 := rec(
        info := "I x L x I loop2",
        forTransposition := false,
        applicable := nt -> not nt.hasTags(),
        apply := (nt, c, cnt) -> let(
            m := nt.params[2], n := nt.params[1]/nt.params[2], j:=Ind(m), fid := fId(n), fbase := fBase(m,j),
            gath := Gath(fTensor(fbase, fid)), scat := Scat(fTensor(fid, fbase)),
            c0 := [ISum(j, m, scat*gath)],
            c1 := When(nt.params[3]=1, [], [I(nt.params[3])]),
            c2 := When(nt.params[4]=1, [], [I(nt.params[4])]),
            Tensor(Concat(c1,c0,c2))
        ),

#D        isApplicable := P -> Length(P[5]) = 0,
#D        rule := (P, C) -> let(m:=P[2], n:=P[1]/P[2], j:=Ind(n), fid := fId(m), fbase := fBase(n,j),
#D                gath := Gath(fTensor(fbase, fid)), scat := Scat(fTensor(fid, fbase)),
#D                C0 := [ISum(j, n, scat*gath)], C1:=When(P[3]=1, [], [I(P[3])]), C2:=When(P[4]=1, [], [I(P[4])]), Tensor(Concat(C1, C0, C2))),
        switch := false
    )
));

###################################################################
NewRulesFor(TICompose, rec(
    TICompose_unroll := rec(
        forTransposition := false,
        applicable := nt -> true,
        children := nt -> [[
            TCompose(
                List([0..nt.params[2]-1], i -> RulesStrengthReduce(SubstBottomUp(Copy(nt.params[3]), nt.params[1], e -> V(i))))
            ).withTags(nt.getTags())
        ]],
        apply := (nt, c, cnt) -> c[1]
    )
));

NewRulesFor(TDR, rec(
    TDR_base := rec(
        forTransposition := false,
        applicable := nt -> true,
        apply := (nt, c, cnt) -> DR(nt.params[1], nt.params[2])
#D        isApplicable := True,
#D        rule := (P, C) -> DR(P[1], P[2])
    )
));

NewRulesFor(TGath, rec(
    TGath_base := rec(
        applicable := True,
        apply := (t, C, nt) -> t.terminate()
    )
));

NewRulesFor(TScat, rec(
    TScat_base := rec(
        applicable := True,
        apply := (t, C, nt) -> t.terminate()
    )
));

NewRulesFor(TConj, rec(
    TConj_tag := rec(
        applicable := True,
        children := t -> [[ t.params[1].withTags(t.getTags()) ]],
        apply := (t, C, nt) -> ConjLR(C[1], t.params[2], t.params[3])
    ),

    TConj_perm := rec(
        applicable := True,

	_cvtPerm := (t,p, use_tl) -> Cond(
	    ObjId(p) = fId,
	        I(p.params[1]),
	    ObjId(p) = L and use_tl,
	        TL(p.params[1], p.params[2], 1, 1).withTags(t.getTags()),
	    # else
		FormatPrm(p)
	),

	# one degree of freedom -- use TL (true) or use FormatPrm(L) (false)
	freedoms := (self, t) >> [[ true, false ]],

        child := (self, t, fr) >> [
	    self._cvtPerm(t, t.params[2], fr[1]),
	    t.params[1].withTags(t.getTags()),
	    self._cvtPerm(t, t.params[3], fr[1])
	],

	apply := (self, t, C, Nonterms) >> C[1]*C[2]*C[3]
    ),

    TConj_cplx := rec(
        applicable := t -> t.params[1] _is TRC,

	_cvtPerm := (t, p) -> Cond(
	    ObjId(p) = fId,
	        I(p.params[1]),
	    ObjId(p) = L,
	        TL(p.params[1], p.params[2], 1, 1).withTags(t.getTags()),
	    # else
		FormatPrm(p)
	),

	freedoms := (self, t) >> [],

        child := (self, t, fr) >> [
	    self._cvtPerm(t, t.params[2]),
	    t.params[1].withTags(List(t.getTags(), t->Cond(t.kind()=spiral.paradigms.vector.AVecReg, spiral.paradigms.vector.AVecRegCx(t.isa.cplx()), t))),
	    self._cvtPerm(t, t.params[3])
	],

	apply := (self, t, C, Nonterms) >> C[1]*C[2]*C[3]
    ),

));

#########################################################################

NewRulesFor(TTensorInd, rec(
#   base cases
#   I x A
    dsA_base := rec(
        info := "IxA base",
        forTransposition := false,
        applicable := nt -> not nt.hasTags() and IsParPar(nt.params),
        children := nt -> [[ nt.params[1], InfoNt(nt.params[2]) ]],
        apply := (nt, c, cnt) -> IDirSum(cnt[2].params[1], c[1])
    ),
#   A x I
    L_dsA_L_base := rec(
        info := "AxI base",
        forTransposition := false,
        applicable := nt -> not nt.hasTags() and IsVecVec(nt.params),
        children := nt -> [[ nt.params[1], InfoNt(nt.params[2]) ]],
        apply := (nt, c, cnt) ->
            L(c[1].dims()[1] * nt.params[2].range, c[1].dims()[1]) *
            IDirSum(cnt[2].params[1], c[1]) *
            L(c[1].dims()[2] * nt.params[2].range, nt.params[2].range)
    ),
#   (I x A)L
    dsA_L_base := rec(
        info := "(IxA)L base",
        forTransposition := false,
        applicable := nt -> not nt.hasTags() and IsParVec(nt.params),
        children := nt -> [[ nt.params[1], InfoNt(nt.params[2]) ]],
        apply := (nt, c, cnt) ->
            IDirSum(cnt[2].params[1], c[1]) *
            L(c[1].dims()[2] * nt.params[2].range, nt.params[2].range),
    ),
#   L(I x A)
    L_dsA_base := rec(
        info := "L(IxA) base",
        forTransposition := false,
        applicable := nt -> not nt.hasTags() and IsVecPar(nt.params),
        children := nt -> [[ nt.params[1], InfoNt(nt.params[2]) ]],
        apply := (nt, c, cnt) ->
            L(c[1].dims()[1] * nt.params[2].range, c[1].dims()[1]) *
            IDirSum(cnt[2].params[1], c[1])
    )
));


text:
--
-- In a parser using this template, the following macro may be redefined:
--
--     $additional_interfaces
--     $ast_class
--
-- B E G I N N I N G   O F   T E M P L A T E   btParserTemplateF
--
%Options programming_language=csharp,margin=4,backtrack
%Options table,error_maps,scopes
%options prefix=TK_
%options action-block=("*.cs", "/.", "./")
%options ParseTable=LPG2.Runtime.ParseTable
%options nt-check

--
-- This template requires that the name of the EOF token be set
-- to EOF_TOKEN to be consistent with LexerTemplateD and LexerTemplateE
--
%EOF
    EOF_TOKEN
%End

%ERROR
    ERROR_TOKEN
%End

%Define

    $Header
    /.
                //
                // Rule $rule_number:  $rule_text
                //
                ./

    $BeginAction
    /.$Header$case $rule_number: {
                   //#line $next_line "$input_file$"./

    $EndAction
    /.            break;
                }./

    $BeginJava
    /.$Header$case $rule_number: {
                    $symbol_declarations
                    //#line $next_line "$input_file$"./

    $EndJava /.$EndAction./

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

    $BadAction
    /.$Header$case $rule_number:
                    throw ("No action specified for rule " + $rule_number);./

    $NullAction
    /.$Header$case $rule_number:
                    setResult(null);
                    break;./

    $BeginActions
    /.
         // Casting object to various generic types
        public void ruleAction(int ruleNumber)
        {
            switch (ruleNumber)
            {./

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

            public void ruleAction$rule_number(int ruleNumber)
            {
                switch (ruleNumber)
                {
                    //#line $next_line "$input_file$"./

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

    $entry_declarations
    /.
        public $ast_class parse$entry_name()
        {
            return parse$entry_name(null, 0);
        }

        public $ast_class parse$entry_name(Monitor monitor)
        {
            return parse$entry_name(monitor, 0);
        }

        public $ast_class parse$entry_name(int error_repair_count)
        {
            return parse$entry_name(null, error_repair_count);
        }

        public $ast_class parse$entry_name(Monitor monitor, int error_repair_count)
        {
            btParser.setMonitor(monitor);

            try
            {
                return ($ast_class) btParser.fuzzyParseEntry($sym_type.$entry_marker, error_repair_count);
            }
            catch (BadParseException e)
            {
                prsStream.reset(e.error_token); // point to error token

                DiagnoseParser diagnoseParser = new DiagnoseParser(prsStream, prsTable);
                diagnoseParser.diagnoseEntry($sym_type.$entry_marker, e.error_token);
            }

            return null;
        }
    ./

    --
    -- Macros that may be needed in a parser using this template
    --
    $additional_interfaces /../
    $ast_class /.$ast_type./
    $super_class /.object./
    $unimplemented_symbols_warning /.false./

    --
    -- Old deprecated macros that should NEVER be used.
    --
    $setSym1 /. // macro setSym1 is deprecated. Use function setResult
                getParser().setSym1./
    $setResult /. // macro setResult is deprecated. Use function setResult
                 getParser().setSym1./
    $getSym /. // macro getSym is deprecated. Use function getRhsSym
              getParser().getSym./
    $getToken /. // macro getToken is deprecated. Use function getRhsTokenIndex
                getParser().getToken./
    $getIToken /. // macro getIToken is deprecated. Use function getRhsIToken
                 prsStream.getIToken./
    $getLeftSpan /. // macro getLeftSpan is deprecated. Use function getLeftSpan
                   getParser().getFirstToken./
    $getRightSpan /. // macro getRightSpan is deprecated. Use function getRightSpan
                    getParser().getLastToken./
%End

%Globals
    /.
    using  LPG2.Runtime;
    using System;
    ./
%End

%Headers
    /.
    public class $action_type : $super_class , RuleAction$additional_interfaces
    {
        private PrsStream prsStream = null;

        private bool unimplementedSymbolsWarning = $unimplemented_symbols_warning;

        private static ParseTable prsTable = new $prs_type();
        public ParseTable getParseTable() { return prsTable; }

        private BacktrackingParser btParser = null;
        public BacktrackingParser getParser() { return btParser; }

        private void setResult(object _object) { btParser.setSym1(_object); }
        public object getRhsSym(int i) { return btParser.getSym(i); }

        public int getRhsTokenIndex(int i) { return btParser.getToken(i); }
        public IToken getRhsIToken(int i) { return prsStream.getIToken(getRhsTokenIndex(i)); }

        public int getRhsFirstTokenIndex(int i) { return btParser.getFirstToken(i); }
        public IToken getRhsFirstIToken(int i) { return prsStream.getIToken(getRhsFirstTokenIndex(i)); }

        public int getRhsLastTokenIndex(int i) { return btParser.getLastToken(i); }
        public IToken getRhsLastIToken(int i) { return prsStream.getIToken(getRhsLastTokenIndex(i)); }

        public int getLeftSpan() { return btParser.getFirstToken(); }
        public IToken getLeftIToken()  { return prsStream.getIToken(getLeftSpan()); }

        public int getRightSpan() { return btParser.getLastToken(); }
        public IToken getRightIToken() { return prsStream.getIToken(getRightSpan()); }

        public int getRhsErrorTokenIndex(int i)
        {
            int index = btParser.getToken(i);
            IToken err = prsStream.getIToken(index);
            return (err is ErrorToken ? index : 0);
        }
        public ErrorToken getRhsErrorIToken(int i)
        {
            int index = btParser.getToken(i);
            IToken err = prsStream.getIToken(index);
            return (ErrorToken) (err is ErrorToken ? err : null);
        }

        public void reset(ILexStream lexStream)
        {
            prsStream = new PrsStream(lexStream);
            btParser.reset(prsStream);

            try
            {
                prsStream.remapTerminalSymbols(orderedTerminalSymbols(), prsTable.getEoftSymbol());
            }
            catch (NullExportedSymbolsException e) {
            }
            catch (NullTerminalSymbolsException e) {
            }
            catch (UnimplementedTerminalsException e)
            {
                if (unimplementedSymbolsWarning) {
                    ArrayListHelper<int> unimplemented_symbols =  new ArrayListHelper<int>( e.getSymbols());
                    Console.Out.WriteLine("The Lexer will not scan the following token(s):");
                    for (int i = 0; i < unimplemented_symbols.Count; i++)
                    {
                        int  id = unimplemented_symbols.get(i);
                        Console.Out.WriteLine("    " + $sym_type.orderedTerminalSymbols[id]);
                    }
                    Console.Out.WriteLine();
                }
            }
            catch (UndefinedEofSymbolException e)
            {
                throw (new UndefinedEofSymbolException
                                    ("The Lexer does not implement the Eof symbol " +
                                     $sym_type.orderedTerminalSymbols[prsTable.getEoftSymbol()]));
            }
        }

        public $action_type()
        {
            try
            {
                btParser = new BacktrackingParser(prsStream, prsTable, (RuleAction) this);
            }
            catch (NotBacktrackParseTableException e)
            {
                throw (new NotBacktrackParseTableException
                                    ("Regenerate $prs_type.cs with -BACKTRACK option"));
            }
            catch (BadParseSymFileException e)
            {
                throw (new BadParseSymFileException("Bad Parser Symbol File -- $sym_type.cs"));
            }
        }

        public $action_type(ILexStream lexStream):this()
        {

            reset(lexStream);
        }

        public int numTokenKinds() { return $sym_type.numTokenKinds; }
        public string[] orderedTerminalSymbols() { return $sym_type.orderedTerminalSymbols; }
        public string getTokenKindName(int kind) { return $sym_type.orderedTerminalSymbols[kind]; }
        public int getEOFTokenKind() { return prsTable.getEoftSymbol(); }
        public IPrsStream getIPrsStream() { return prsStream; }

        /**
         * @deprecated replaced by {@link #getIPrsStream()}
         *
         */
        public PrsStream getPrsStream() { return prsStream; }

        /**
         * @deprecated replaced by {@link #getIPrsStream()}
         *
         */
        public PrsStream getParseStream() { return prsStream; }

        public $ast_class parser()
        {
            return parser(null, 0);
        }

        public $ast_class parser(Monitor monitor)
        {
            return parser(monitor, 0);
        }

        public $ast_class parser(int error_repair_count)
        {
            return parser(null, error_repair_count);
        }

        public $ast_class parser(Monitor monitor, int error_repair_count)
        {
            btParser.setMonitor(monitor);

            try
            {
                return ($ast_class) btParser.fuzzyParse(error_repair_count);
            }
            catch (BadParseException e)
            {
                prsStream.reset(e.error_token); // point to error token

                DiagnoseParser diagnoseParser = new DiagnoseParser(prsStream, prsTable);
                diagnoseParser.diagnose(e.error_token);
            }

            return null;
        }

        //
        // Additional entry points, if any
        //
        $entry_declarations
    ./

%End

%Rules
    /.$BeginActions./
%End

%Trailers
    /.
        $EndActions
    }
    ./
%End

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


text:
#
# Read("~/Workspace/groupsSB/epi/A2/handle.gi");
#

type:="A";
rank:=2;
nr_pos_roots:=3;

Read("~/Workspace/groupsSB/epi/group.gi");

#sList(rels,r->Value(r,[vars[22]],[vars[2]+vars[12]]);

#
# c_1, c_2 not 0 , so they can be choosen c_1=c_2=1 => c_3=c_1/2
#
handleReg:=function()
	# a test that the coeffs are right
	#a:=root_group(1,xvars[1]);
	#b:=root_group(2,xvars[2]);
	#c:=root_group(3,xvars[3]);
	#test1:=evaluate_U(a*b*c,[[[xvars[1]],[One(APR)]],[[xvars[2]],[One(APR)]],[[xvars[3]],[One(APR)/2]]]);
	#test2:=evaluate_U(a*b*c,[[[xvars[1]],[One(APR)*2]],[[xvars[2]],[One(APR)*2]],[[xvars[3]],[One(APR)*2]]]);
	#test2-test1^2;

	local ua,ub,uab,u,e,nn,vals,rels,i,v;
	#ua:=root_group(1,2);
	#ub:=root_group(2,2);
	#uab:=root_group(3,2);
	ua:=root_group(1,1);
	ub:=root_group(2,1);
	uab:=root_group(3,1/2);
	u:=ua*ub*uab;
	e:=u-id_mat;

	nn:=List([1..2*nr_pos_roots+rank],i->xvars[i]);

	vals:=[];
	Append(vals,[[[xvars[6]],[Zero(APR)]]]);
	Append(vals,[[[xvars[4]],[Zero(APR)]]]);
	Append(vals,[[[xvars[5]],[Zero(APR)]]]);
	Append(vals,[[[xvars[7]],[2*xvars[8]]]]);
	Append(vals,[[[xvars[8]],[Zero(APR)]]]);
	Append(vals,[[[xvars[2]],[xvars[1]]]]);
	##
	#Append(vals,[[[xvars[1]],[Zero(APR)]]]);
	#Append(vals,[[[xvars[3]],[Zero(APR)]]]);
	#Append(vals,[[[xvars[7]],[Zero(APR)]]]);

	#
	# atentie la ordinea transpunerilor
	#
	rels:=nn*TransposedMat(u)-nn;
	for i in [1..Length(rels)] do
		for v in vals do
			Print(rels[i],"\n");
			rels[i]:=One(APR)*Value(rels[i],v[1],v[2]);
			nn[i]:=One(APR)*Value(nn[i],v[1],v[2]);
		od;
	od;
	return [u,nn,rels];
end;

#
# c_2=0
#
handle2:=function()
	local ua,ub,uab,u,e,nn,vals,rels,i,v;
	ua:=root_group(1,1);
	uab:=root_group(3,1);
	u:=ua*uab;
	e:=u-id_mat;

	nn:=List([1..2*nr_pos_roots+rank],i->xvars[i]);

	vals:=[];
	Append(vals,[[[xvars[4]],[Zero(APR)]]]);
	Append(vals,[[[xvars[6]],[Zero(APR)]]]);
	Append(vals,[[[xvars[2]],[-xvars[7]-xvars[8]]]]);
	Append(vals,[[[xvars[5]],[-2*xvars[7]+xvars[8]]]]);
	##
	Append(vals,[[[xvars[1]],[Zero(APR)]]]);
	Append(vals,[[[xvars[3]],[Zero(APR)]]]);
	Append(vals,[[[xvars[7]],[Zero(APR)]]]);

	#
	# atentie la ordinea transpunerilor
	#
	rels:=nn*TransposedMat(u)-nn;
	for i in [1..Length(rels)] do
		for v in vals do
			Print(rels[i],"\n");
			rels[i]:=One(APR)*Value(rels[i],v[1],v[2]);
			nn[i]:=One(APR)*Value(nn[i],v[1],v[2]);
		od;
	od;
	return [u,nn,rels];
end;

max_stars_repo_name: iuliansimion/groupsSB

max_issues_repo_name: iuliansimion/groupsSB

max_forks_repo_name: iuliansimion/groupsSB


text:
#############################################################################
####
##
#W  ace.gi                     ACE Package                   Alexander Hulpke
#W                                                                Greg Gamble
##
##  `Head' file for the GAP interface to the ACE (Advanced Coset Enumerator),
##  by George Havas and Colin Ramsay.  The original interface was written  by
##  Alexander Hulpke and extensively modified by Greg Gamble.
##
#Y  Copyright (C) 2000  Centre for Discrete Mathematics and Computing
#Y                      Department of Information Technology & Electrical Eng.
#Y                      University of Queensland, Australia.
##

#############################################################################
####
##
#V  ACEData . . . . . . . record used by various functions of the ACE package
##
##  The fields of ACEData are:
##
##    "binary"  . . the path of the ACE binary
##    "tmpdir"  . . the path of the temporary directory for ACE i/o files
##    "ni"  . . . . record for a non-interactive process
##    "io"  . . . . list of data records for ACEStart IO Streams
##    "infile"  . . the path of the ACE input file
##    "outfile" . . the path of the ACE output file
##    "version" . . the version of the current ACE binary
##
InstallValue( ACEData,
  rec( binary := ExternalFilename(DirectoriesPackagePrograms("ace"), "ace"),
       tmpdir := DirectoryTemporary(),
       ni     := rec(),
       io     := [] # Initially no ACEStart IO Streams
       )
);
ACEData.infile  := Filename(ACEData.tmpdir, "in");
ACEData.outfile := Filename(ACEData.tmpdir, "out");

PrintTo(ACEData.infile, "\n");
# Fire up ACE with a null input (ACEData.infile contains only a "\n")
# ... to generate a banner (which has ACE's current version)
Exec(Concatenation(ACEData.binary, "<", ACEData.infile, ">", ACEData.outfile));
ACEData.version := StringFile( ACEData.outfile );
ACEData.scratch := PositionSublist(ACEData.version, "ACE") + 4;
ACEData.version := ACEData.version{[ACEData.scratch ..
                                    Position(ACEData.version, ' ',
                                             ACEData.scratch) - 1]};
Unbind(ACEData.scratch); # We don't need ACEData.scratch, anymore.

#############################################################################
##
#I  InfoClass
##
# Set the default level of InfoACE
SetInfoLevel(InfoACE, 1);

#############################################################################
####
##
#V  ACEIgnoreUnknownDefault . . . . . . . . . . . .  the default value of the
##  . . . . . . . . . . . . . . . . . . . . . . . . `aceignoreunknown' option
##
ACEIgnoreUnknownDefault := true;

#############################################################################
####
##  Ensure no zombie ACE processes from interactive (ACEStart)  sessions  are
##  . . . . .  . . . . . . . . . . . .  left lying around when user quits GAP
##
InstallAtExit( ACEQuitAll );

#E  ace.gi . . . . . . . . . . . . . . . . . . . . . . . . . . . .  ends here

max_stars_repo_name: isuruf/ace

max_issues_repo_name: isuruf/ace

max_forks_repo_name: isuruf/ace


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

Class(MultiPtrCodegenMixin, rec(
    GathPtr := meth(self, o, y, x, opts)
        local i, func, rfunc, ix;
        i := Ind(); func := o.func.lambda();

        o.ptr.setType();
        if IsBound(o.func.rlambda) then
            rfunc := o.func.rlambda();
            ix := var.fresh_t("ix", TInt);
            return decl(ix, chain(
                    assign(ix, func.at(0)),
                    assign(nth(y, 0), nth(o.ptr, ix)),
                    loop(i, o.func.domain()-1,
                        chain(assign(ix, rfunc.at(ix)),
                            assign(nth(y, i+1), nth(o.ptr, ix))))));
        else
            return loop(i, o.func.domain(), assign(nth(y,i), nth(o.ptr, func.at(i))));
        fi;
    end,

    ScatPtr := meth(self, o, y, x, opts)
        local i, func, rfunc, ix;
        i := Ind(); func := o.func.lambda();
        o.ptr.setType();

        if IsBound(o.func.rlambda) then
            rfunc := o.func.rlambda();
            ix := var.fresh_t("ix", TInt);
            return decl(ix, chain(
                    assign(ix, func.at(0)),
                    assign(nth(o.ptr, ix), nth(x, 0)),
                    loop(i, o.func.domain()-1,
                        chain(assign(ix, rfunc.at(ix)),
                            assign(nth(o.ptr, ix), nth(x, i+1))))));
        else
            return loop(i, o.func.domain(), assign(nth(o.ptr, func.at(i)), nth(x, i)));
        fi;
    end
));

max_stars_repo_name: franzfranchetti/spiral-package-fftx

max_issues_repo_name: franzfranchetti/spiral-package-fftx

max_forks_repo_name: franzfranchetti/spiral-package-fftx


text:
RequirePackage("grape");

# bilinear forms graph: vozlišča so matrike,
# dve vozlišči sta sosednji, ko je njuna razlika ranga 1
GB := AdjFunGraph(Elements(GF(2)^[2, 3]),
    function(x,y) return RankMat(x-y) = 1; end);;

# OA(8, 3): vozlišča so vnosi v latinskem kvadratu (Z_8)^2,
# dve vozlišči sta sosednji, če sta v isti vrstici, v istem stolpcu,
# ali pa imata isto vrednost
GZ8 := AdjFunGraph(Elements(ZmodnZ(8)^2), function(x,y) return x <> y and
        (x[1] = y[1] or x[2] = y[2] or x[1]+x[2] = y[1]+y[2]); end);;

# OA(8, 3) za latinski kvadrat (Z_2 × Z_4)^2
A := Elements(AbelianGroup([2, 4]));;
A2 := Cartesian(A, A);;
GZ24 := AdjFunGraph(A2, function(x,y) return x <> y and
        (x[1] = y[1] or x[2] = y[2] or x[1]*x[2] = y[1]*y[2]); end);;

# OA(8, 3) za latinski kvadrat (Dih(8))^2
D := Elements(DihedralGroup(8));;
D2 := Cartesian(D, D);;
GD8 := AdjFunGraph(D2, function(x,y) return x <> y and
        (x[1] = y[1] or x[2] = y[2] or x[1]*x[2] = y[1]*y[2]); end);;

# OA(8, 3) za latinski kvadrat (Q_8)^2
Q := Elements(SmallGroup(8, 4));;
Q2 := Cartesian(Q, Q);;
GQ8 := AdjFunGraph(Q2, function(x,y) return x <> y and
        (x[1] = y[1] or x[2] = y[2] or x[1]*x[2] = y[1]*y[2]); end);;

# iz [21, 6; 8:21, 12:42]_2 kode
M := [
    [1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0],
    [1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0],
    [1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0],
    [1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0],
    [1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 0],
    [1, 0, 1, 0, 1, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1]
]*Z(2);

# vozlišča so kodne besede,
# dve vozlišči sta sosednji, ko je razdalja med njima minimalna
GC := AdjFunGraph(Elements(Subspace(GF(2)^21, M)),
    function(x, y) return WeightVecFFE(x+y) = 8; end);;
lGC := DistanceSetInduced(GC, 1, 1);;

ext: gap

max_stars_repo_name: jaanos/gap-graphs

max_issues_repo_name: jaanos/gap-graphs

max_forks_repo_name: jaanos/gap-graphs


text:
f := function(x,y)
    local z, helper;
    z := 1 - y;
    helper := function(bla)
        return x + y*z - bla;
    end;
    return helper;
end;

x -> x + 1;
{x,y} -> x + y + 1;

Print("This is a string\n");

for n in [1..100] do
  for i in [1..NrSmallGroups(n)] do
    G := SmallGroup(n,i);
    Print("SmallGroup(", n, ",", i, ") is abelian? ", IsAbelian(G), "\n");
    if IsPerfectGroup(G) = true then
      break;
    fi;
  od;
od;

BindGlobal("foo", 1);
BIND_GLOBAL("bar", 2);

r := rec(
    a := [1,2,3],
    b := true,
    c := fail,
   );

InstallImmediateMethod( IsFinitelyGeneratedGroup,
    IsGroup and HasGeneratorsOfGroup,
    function( G )
    if IsFinite( GeneratorsOfGroup( G ) ) then
      return true;
    fi;
    TryNextMethod();
    end );

max_stars_repo_name: fingolfin/tree-sitter-gap

max_issues_repo_name: fingolfin/tree-sitter-gap

max_forks_repo_name: fingolfin/tree-sitter-gap


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

TestFailExit := function()
    Print("TEST FAILED\n");
    Exit(-1);
end;

TestSkipExit := function()
    Print("Skipping test\n");
    Exit(86);
end;

#F  GetBasicProfilerTestFName (<pf>)
#F      Get the filename associated with Basic Profiler Test
#F      if <pf> return Success file name
#F      otherwise return Failure file name

GetBasicProfilerTestFName := function(pf)
    local path, sep;
    sep  := Conf("path_sep");
    path := Conf("spiral_dir");
    path := Concat(path, sep, "build");
    if pf then
        path := Concat(path, sep, "PROFILER_RUN_SUCCESS");
    else
        path := Concat(path, sep, "PROFILER_RUN_FAILED");
    fi;
    return path;
end;

#F  ClearBasicProfilerTestResults () -- remove results of prior run(s) of basic profiler test

ClearBasicProfilerTestResults := function()
    SysRemove(GetBasicProfilerTestFName(true));
    SysRemove(GetBasicProfilerTestFName(false));
    return;
end;

#F  MarkBasicProfilerTest (<pf>)
#F      Mark (create the filename) associated with Basic Profiler Test
#F      if <pf> ==> test successed, use success file name
#F      otherwise ==> test failed, use Failure file name

MarkBasicProfilerTest := function(pf)
    local path;
    path := GetBasicProfilerTestFName(pf);
    PrintTo(path, "");
    return;
end;

#F  CheckBasicProfilerTest () -- Return True if basic profiler test passed, otherwise, False

CheckBasicProfilerTest := function()
    local res, file;
    file := GetBasicProfilerTestFName(true);
    res  := CheckFileExists(file, "");

    if SysVerbose() > 0 then
        Print("Marker file: ", file);
        if res then PrintLine(" Found, return true"); else PrintLine(" NOT Found"); fi;
    fi;
    if res then return res; fi;

    file := GetBasicProfilerTestFName(false);
    res := CheckFileExists(file, "");
    if SysVerbose() > 0 then
        Print("Marker file: ", file);
        if res then PrintLine(" Found, return false"); else PrintLine(" NOT Found"); fi;
    fi;
    return false;
end;


text:
# Examples/Tests for ANATPH

# Jack Button's Group
# -------------------
#
# With presentation
#
# Alan Logan says this is hyperbolic, but noone wants to publish
# this result alone.
#
# Tester now proves this to be hyperbolic, after corrected
# relation. Bug reported by email by Chris Chalk <chalk235@gmail.com>
#
InstallGlobalFunction("JackButtonGroup",
function()
    local pg;

    pg := PregroupOfFreeGroup(3);
    SetPregroupElementNames(pg, "1aAbBtT");
    return NewPregroupPresentation(pg, [ pg_word( pg, [7,2,6,5,3])
                                       , pg_word( pg, [7,4,6,3,5]) ]);
end);

# Triangle Groups
# ---------------
# As demonstrated in Proposition 9.4 of anatph, for l > 2 (and hence n > 3 if the
# group is supposed to be hyperbolic) there are no instantiable green places, and if
# l = 2 then there are exactly two instantiable green places.
InstallGlobalFunction("TriangleGroup",
function(l,m,n)
    local pg;
    pg := PregroupOfFreeProduct( CyclicGroup( IsPermGroup, l)
                               , CyclicGroup( IsPermGroup, m) );

    # This is a bit icky, we can't tell which element of pg is the one of the
    # second cyclic group
    return NewPregroupPresentation(pg, [ pg_word( pg, Repeat(n, [2, l + 1]))]);
end);

# Example from Theorem 9.5, triangle-like: quotient of 2-3-m triangle group
# our choice of pregroup is the free product of cyclic groups of order 2 and 3
#T see what happens if we present this group as
#T <x,y,z | x^2, y^3, z^m, (zxY)^n > ?
#T we push a parameter into the pregroup
InstallGlobalFunction(TriangleCommutatorQuotient,
function(m,n)
    local pg;
    pg := PregroupOfFreeProduct( CyclicGroup(IsPermGroup, 2)
                               , CyclicGroup(IsPermGroup, 3) );
    # Do this to have slightly nicer display. Maybe we need to give the user
    # a way to label generators
    pg!.enams := "1xyY";
    return NewPregroupPresentation(pg,
                                   [ pg_word(pg, Repeat(m, [2,3])),
                                     pg_word(pg, Repeat(n, [2,3,2,4]))
                                   ]);
end);

InstallGlobalFunction(RandomTriangleQuotient,
function(p,q,r,len)
    local pg;
    pg := PregroupOfFreeProduct( CyclicGroup(IsPermGroup, p)
                               , CyclicGroup(IsPermGroup, q) );
    # Do this to have slightly nicer display. Maybe we need to give the user
    # a way to label generators
    pg!.enams := "1xyY";
    return NewPregroupPresentation(pg,
                                   [ pg_word(pg, Repeat(r, [2,3])),
                                     RandomPregroupWord(pg, len)
                                   ]);
end);

InstallGlobalFunction(RandomPregroupWord,
function(pg, len)
    local i, lett, rel;

    rel := [];

    rel[1] := Random([2..Size(pg)]);
    for i in [2..len-1] do
        rel[i] := Random(Difference( [2..Size(pg)]
                                   , [__ID(PregroupInverse(pg[rel[i-1]]))]));
    od;
    rel[len] := Random(Difference([2..Size(pg)]
                                 , [ __ID(PregroupInverse(pg[rel[len - 1]]))
                                   , __ID(PregroupInverse(pg[rel[1]]))
                                   ] ));
    return pg_word(pg, rel);
end);

# Given a pregroup make a random presentation with nrel relators
# of length lrel
InstallGlobalFunction(RandomPregroupPresentation,
function(pg, nrel, lrel)
    local rels;

    rels := List([1..nrel], i -> RandomPregroupWord(pg, lrel));
    return NewPregroupPresentation(pg, rels);
end);

BindGlobal("CreateRandomExample",
function(path, pg, nrel, lrel)
    local pgp;

    if not IsDirectoryPath(path) then
        Error("path does not exist or is not a directory");
        return;
    fi;
    pgp := RandomPregroupPresentation(pg, nrel, lrel);

    # Write the pregroup presentation to file
    PregroupPresentationToFile(Concatenation(path, "/presentation-gap"), pgp);
    PregroupPresentationToSimpleFile(Concatenation(path, "/presentation-simple"), pgp);

    if WALRUS_kbmag_available then
    # Also create a KBMAG input file, if RSymTest succeeds on this presentation,
    # we can try computing an automatic structure and run gpgeowa on the result
    # to check whether the group is hyperbolic.
    WriteRWS(KBMAGRewritingSystem(PregroupPresentationToFpGroup(pgp)),
             Concatenation(path, "/presentation-kbmag"));
    fi;
end);

# Create a series of examples
BindGlobal("CreateRandomSeries",
function(pg, nrels, lrels, nexs, prf, basepath)
    local path, path2, i, pgp, lrel, nrel;

    if not IsDirectoryPath(basepath) then
        Error("path does not exist or is not a directory");
        return;
    fi;
    path := basepath;

    for lrel in lrels do
        CreateDir(Concatenation(path, "/", String(lrel)));
        for nrel in nrels do
            CreateDir(Concatenation(path, "/", String(lrel), "/", String(nrel)));
            for i in [1..nexs] do
                path2 := Concatenation(path, "/", String(lrel), "/", String(nrel), "/", String(i));
                prf("creating example nr ", i, " in ", path2, " \c");

                CreateDir(path2);
                CreateRandomExample(path2, pg, nrel, lrel);
                prf("\n");
            od;
        od;
    od;
end);

BindGlobal("CreateRandomSeriesOverSmallPregroups",
function(basepath)
    local i, n, sizes, path;

    if not IsDirectoryPath(basepath) then
        Error("path does not exist or is not a directory");
        return;
    fi;

    # at the moment we only have pregroups
    # of size 6
    sizes := [6];
    for n in sizes do
        CreateDir(Concatenation(basepath, "/", String(n)));
        for i in [1..NrSmallPregroups(n)] do
            path := Concatenation(basepath, "/", String(n), "/", String(i));
            CreateDir(path);
            CreateRandomSeries( SmallPregroup(n, i)
                              , [1,2,4,8,10]
                              , [5,10,15,20,25,30,35,40,45,50]
                              , 20
                              , Print
                              , path);
        od;
    od;
end);

BindGlobal("CreateRandomSeriesOverFreeGroup",
function(basepath)
    local i, n, sizes, path;

    if not IsDirectoryPath(basepath) then
        Error("path does not exist or is not a directory");
        return;
    fi;

    # ranks of free group
    sizes := [2,4,8,16,32,64];
    for n in sizes do
        path := Concatenation(basepath, "/", String(n));
        CreateDir(path);
        CreateRandomSeries( PregroupOfFreeGroup(n)
                          , [1,2,4,8,10]
                          , [5,10,15,20,25,30,35,40,45,50]
                          , 20
                          , Print
                          , path);
    od;
end);

BindGlobal("CreateRandomSeriesOverTrianglePregroup",
function(eps, p, q, nrel, lrel, nexs, prf, path)
    local pg;
    pg := PregroupOfFreeProduct( CyclicGroup( IsPermGroup, p)
                                , CyclicGroup( IsPermGroup, q) );
    pg!.enams := "1xyY";
    CreateRandomSeries( pg
                      , [1]
                      , [5,10,15,20,25,30,35,40,45,50]
                      , 20
                      , Print
                      , path );
end);

BindGlobal("CreateRandomSeriesOverTriangleGroup",
          function(eps, p, q, nrel, lrel, nexs, prf, path)
              local pg;
              pg := PregroupOfFreeProduct( CyclicGroup( IsPermGroup, p)
                                         , CyclicGroup( IsPermGroup, q) );
              pg!.enams := "1xyY";
              CreateRandomSeries( pg
                                , [1]
                                , [5,10,15,20,25,30,35,40,45,50]
                                , 20
                                , Print
                                , path );
          end);

BindGlobal("BenchmarkRandomPresentation",
function(pg, eps, nrel, lrel, nexs, prf, path)
    local i, pgp, start, stop, estart, estop, n, nfail, res, runt, fid, stream;

    n := 0;
    runt := [];

    start := NanosecondsSinceEpoch();
    for i in [1..nexs] do
        prf("creating example nr ", i, ", \c");
        pgp := RandomPregroupPresentation(pg, nrel, lrel);
        prf("starting RSymTest, \c");
        estart := NanosecondsSinceEpoch();
        res := RSymTest(pgp, eps);
        estop := NanosecondsSinceEpoch();
        if res = true then
            prf("succeeded after \c");
        else
            prf("failed \c");
        fi;
#        LogPregroupPresentation(path, pgp, res);
        prf(Float((estop - estart) / 1000000000), " seconds\n");
        Add(runt, [pgp, res, estop - estart]);
    od;
    stop := NanosecondsSinceEpoch();
    return rec( runtime := (stop - start) / 1000000000,
                samples := runt );
end);

BindGlobal("RandomPregroupFromSmallGroups",
function()
    local n, i, g1, g2;

    n := Random([1..64]);
    i := Random([1..NrSmallGroups(n)]);
    g1 := SmallGroup(n,i);

    n := Random([1..64]);
    i := Random([1..NrSmallGroups(n)]);
    g2 := SmallGroup(n,i);

    return PregroupOfFreeProduct(g1, g2);
end);

BindGlobal("BenchmarkRandom_TriPregroup",
function(eps, nrel, lrel, nexs, prf, path)
    local pg;
    pg := PregroupOfFreeProduct( CyclicGroup( IsPermGroup, 2)
                               , CyclicGroup( IsPermGroup, 3) );
    pg!.enams := "1xyY";
    return BenchmarkRandomPresentation(pg, eps, nrel, lrel, nexs, prf, path);
end);

BindGlobal("BenchmarkRandom_FreeGroupPregroup",
function(eps, ngen, nrel, lrel, nexs, prf, path)
    local pg;
    pg := PregroupOfFreeGroup(ngen);
    return BenchmarkRandomPresentation(pg, eps, nrel, lrel, nexs, prf, path);
end);

BindGlobal("BenchmarkRandom_OverSmallPregroup",
function(eps, nrel, lrel, nexs, prf, path)
    local pg, res;

    res := [];
    # FIXME: We only have pregroups of size 6 at the moment
    for pg in ANATPH_small_pregroups[6] do
        pg := PregroupByTable([1,'a','b','c','d', 'e'], pg);
        Add(res, BenchmarkRandomPresentation(pg, eps, nrel, lrel, nexs, prf, path));
    od;
    return res;
end);

BindGlobal("_VARIANCE",
function(v)
    local q, avg;
    if Length(v) = 1 then
        q := 1;
    else
        q := Length(v) - 1;
    fi;
    avg := Average(v);
    return Sum(List(v, x -> (x - avg)^2)) / (q);
end);

BindGlobal("AnalyseBenchmarkResult",
function(res)
    local r, succed, failed;

    r := rec();

    r.n := Length(res.samples);
    r.rt_avg := Average(List(res.samples, x->x[3]));
    r.rt_std := Sqrt(Float(_VARIANCE(List(res.samples, x->x[3]))));

    succed := Filtered(res.samples, x -> not IsList(x[2]));
    r.succed_n := Length(succed);
    if r.succed_n > 0 then
        r.succed_rt_avg := Average(List(succed, x->x[3]));
        r.succed_rt_std := Sqrt(Float(_VARIANCE(List(succed, x->x[3]))));
    fi;

    failed := Filtered(res.samples, x -> IsList(x[2]));
    r.failed_n := Length(failed);
    if r.failed_n > 0 then
        r.failed_rt_avg := Average(List(failed, x->x[3]));
        r.failed_rt_std := Sqrt(Float(_VARIANCE(List(failed, x->x[3]))));
    fi;

    return r;
end);

BindGlobal("StringBenchResult",
function(r)
    local res;

    res := STRINGIFY("n: ", r.n, ", s: ", r.succed_n, " f: ", r.failed_n, " (", Float(r.succed_n) / r.n, ")\n",
                     " (", Float(r.rt_avg), ",", r.rt_std, ")\n");
    Append(res, STRINGIFY(" success: ", r.succed_n));
    if r.succed_n > 0 then
        Append(res, STRINGIFY(" ", Float(r.succed_rt_avg), ",", r.succed_rt_std, "\n"));
    fi;
    Append(res, STRINGIFY(" failure: ", r.failed_n));
    if r.failed_n > 0 then
        Append(res, STRINGIFY(" ", Float(r.failed_rt_avg), ",", r.failed_rt_std, "\n"));
    fi;
    if res[Length(res)] <> '\n' then
        Add(res, '\n');
    fi;
    return res;
end);

BindGlobal("StandardBenchmarks",
function(nsamples, prn)
    local d, ngens, rlen, nrels, r, path;

    path := DirectoryTemporary();
    prn("running standard benchmarks for anatph, ", nsamples, " each\n");

    prn("= = = = = = =\n");
    prn(" Quotients of free group pregroup\n");
    prn("= = = = = = =\n");
    for ngens in [2,4,8,16] do
        for nrels in [1,2,3,4] do
            for rlen in [10,20,30,40,50] do
                prn("running "
                   , ngens, " generators, "
                   , nrels, " relators "
                   , "of length ", rlen, "\n");
                r := BenchmarkRandom_FreeGroupPregroup(
                                                    1/12
                                                  , ngens
                                                  , nrels
                                                  , rlen
                                                  , nsamples
                                                  , prn
                                                  , path
                     );
                prn(StringBenchResult(AnalyseBenchmarkResult(r)));
            od;
        od;
    od;

    prn("= = = = = = =\n");
    prn(" Quotients of 2-3-triangle group pregroup\n");
    prn("= = = = = = =\n");
    for nrels in [1,2,3,4] do
        for rlen in [10,20,30,40,50] do
            prn("running ", nrels, "relators of length ", rlen, "\n");
            r := BenchmarkRandom_TriPregroup(
                                              1/12
                                            , nrels
                                            , rlen
                                            , nsamples
                                            , prn
                                            , path
                     );
            prn(StringBenchResult(AnalyseBenchmarkResult(r)));
        od;
    od;

end);

BenchmarkSinglePres := function(eps, pg, nrel, len)
    local r,t,res;
    r := RandomPregroupPresentation(pg, nrel, len);
    t := NanosecondsSinceEpoch();
    res := RSymTest(r, eps);
    t := NanosecondsSinceEpoch() - t;
    return [r, res, t / 1000000000.];
end;

if IsBound(OutputAnnotatedCodeCoverageFiles) then
    ProfileSinglePresentation := function(eps, pg, nrel, len)
    local t, dir, fn;

    dir := DirectoryTemporary();
    fn := Filename(dir, "anatph.gz");

    Print("Profiling presentation to ", fn, "\n");
    ProfileLineByLine(fn);
    t := BenchmarkSinglePres(eps, pg, nrel, len);
    UnprofileLineByLine();

    Print("Writing annotated code coverage to ", Filename(dir, ""), "\n");
    OutputAnnotatedCodeCoverageFiles(fn, Filename(dir, ""));

    return t;
end;

else
    Print("profiling package is not available, disabling ProfileSinglePresentation\n");
fi;

# Benchmark the RSym tester with a presentation given density
# and length
BenchmarkGromovDensity := function(eps, ngens, length, density)
    local pg, nrels;

    pg := PregroupOfFreeGroup(ngens);
    nrels := Int( (2 * ngens - 1)^Float(density * length ) );
    Print("number of relators: ", nrels, "\n");
    return BenchmarkSinglePres(eps
                              , pg
                              , nrels
                              , length );
end;

max_stars_repo_name: RussWoodroofe/walrus

max_issues_repo_name: RussWoodroofe/walrus

max_forks_repo_name: RussWoodroofe/walrus


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

Class(MD_PRDFT, TaggedNonTerminal, rec(
    abbrevs := [ (dims)     -> Checked(IsList(dims), ForAll(dims, IsPosIntSym), [dims, 1]),
                 (dims,rot) -> Checked(IsList(dims), ForAll(dims, IsPosIntSym), IsIntSym(rot), [dims, rot])
    ],

    terminate := self >> let(
        n := self.params[1], rot := self.params[2],
        n0 := DropLast(n, 1),
        lst := Last(n),
        res := RC(Tensor(MDDFT(n0), I(Rows(PRDFT(lst))/2))) * Tensor(I(Product(n0)), PRDFT(lst, rot)),
        When(self.transposed, res.transpose(), res)),

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

    isReal := True,

    normalizedArithCost := self >> let(n := Product(self.params[1]),
       IntDouble(2.5 * n * d_log(n) / d_log(2))),

    hashAs := self >> let(t:=ObjId(self)(self.params[1], 1).withTags(self.getTags()),
        When(self.transposed, t.transpose(), t)),

    dims := self >> let(
        n := self.params[1],
        n0 := DropLast(n, 1),
        lst := Last(n),
        d := [ Product(n0) * 2*(idiv(lst,2)+1), Product(n)],
        When(self.transposed, Reversed(d), d)),

    omega := (N,k,r,c) -> E(N)^(k*r*c),
));

NewRulesFor(MD_PRDFT, rec(
    MD_PRDFT_Base := rec(
        applicable := t -> true,
        freedoms := t -> [],
        child := (t, fr) -> let(
            n := t.params[1], rot := t.params[2],
            n0 := DropLast(n, 1),
            lst := Last(n),
            RC(
                GT(MDDFT(n0), GTVec, GTVec, [Rows(PRDFT(lst))/2]).withTags(t.getTags())
            ) *
            GT(PRDFT(lst, rot), GTPar, GTPar, [Product(n0)]).withTags(t.getTags())
        ),
        apply := (self, t, C, Nonterms) >> C[1]
    )
));

# tst := vec -> let(n := Length(vec), dim := Sqrt(n), cols := 2*(Int(dim/2)+1),
#     out := TransposedMat(MatSPL(MD_PRDFT([dim,dim], 1)) * TransposedMat([vec]))[1],
#     outmat := List([0..dim-1], x->out{[1 + x*cols .. cols + x * cols]}),
#     PrintMat(outmat),
#     outmat);

# ctst := vec -> let(n := Length(vec), dim := Sqrt(n),
#     out := TransposedMat(MatSPL(MDDFT([dim,dim], 1)) * TransposedMat([vec]))[1],
#     outmat := List([0..dim-1], x->out{[1 + x*dim .. dim + x * dim]}),
#     PrintMat((outmat)),
#     outmat);

# CX = X * RE
# X^-1 = RE * CX^-1


text:
#
# francy: Interactive Discrete Mathematics in GAP
#

#############################################################################
##
#M  Canvas( <title>, <options> ) . . . . . a new graphic canvas
##
InstallMethod(Canvas,
  "a title string, a default configurations record",
  true,
  [IsString,
   IsCanvasDefaults],
  0,
function(title, options)
  return MergeObjects(Objectify(CanvasObjectType, rec(
    id       := GenerateID(),
    menus    := rec(),
    graph    := rec(),
    chart    := rec(),
    messages := rec(),
    title    := title
  )), options);
end);

InstallOtherMethod(Canvas,
  "a title string",
  true,
  [IsString],
  0,
function(title)
  return Canvas(title, CanvasDefaults);
end);

#############################################################################
##
#M  Add( <canvas>, <francy object> ) . . . . . add objects to canvas
##
InstallOtherMethod(Add,
  "a canvas, a graph",
  true,
  [IsCanvas,
   IsFrancyGraph],
  0,
function(canvas, graph)
  canvas!.graph := graph;
  # unbind the chart, only graph should exist!
  Unbind(canvas!.chart);
  return canvas;
end);

InstallOtherMethod(Add,
  "a canvas, a chart",
  true,
  [IsCanvas,
   IsChart],
  0,
function(canvas, chart)
  canvas!.chart := chart;
  # unbind the graph, only chart should exist!
  Unbind(canvas!.graph);
  return canvas;
end);

InstallOtherMethod(Add,
  "a canvas, a menu",
  true,
  [IsCanvas,
   IsMenu],
  0,
function(canvas, menu)
  canvas!.menus!.(menu!.id) := menu;
  return canvas;
end);

InstallOtherMethod(Add,
  "a canvas, a message",
  true,
  [IsCanvas,
   IsFrancyMessage],
  0,
function(canvas, message)
  canvas!.messages!.(message!.id) := message;
  return canvas;
end);

InstallOtherMethod(Add,
  "a canvas, a list of francy objects",
  true,
  [IsCanvas,
   IsList],
  0,
function(canvas, objects)
  local object;
  for object in objects do
    Add(canvas, object);
  od;
  return canvas;
end);

#############################################################################
##
#M  Remove( <canvas>, <francy object> ) . . . . . remove object from canvas
##
InstallOtherMethod(Remove,
  "a canvas, a graph",
  true,
  [IsCanvas,
   IsFrancyGraph],
  0,
function(canvas, graph)
  Unbind(canvas!.graph);
  canvas!.graph := rec();
  canvas!.chart := rec();
  return canvas;
end);

InstallOtherMethod(Remove,
  "a canvas, a chart",
  true,
  [IsCanvas,
   IsChart],
  0,
function(canvas, chart)
  Unbind(canvas!.chart);
  canvas!.graph := rec();
  canvas!.chart := rec();
  return canvas;
end);

InstallOtherMethod(Remove,
  "a canvas, a menu",
  true,
  [IsCanvas,
   IsMenu],
  0,
function(canvas, menu)
  Unbind(canvas!.menus!.(menu!.id));
  return canvas;
end);

InstallOtherMethod(Remove,
  "a canvas, a message",
  true,
  [IsCanvas,
   IsFrancyMessage],
  0,
function(canvas, message)
  Unbind(canvas!.messages!.(message!.id));
  return canvas;
end);

InstallOtherMethod(Remove,
  "a canvas, a list of francy objects",
  true,
  [IsCanvas,
   IsList],
  0,
function(canvas, objects)
  local object;
  for object in objects do
    Remove(canvas, object);
  od;
  return canvas;
end);

#############################################################################
##
#M  Draw( ) . . . . .
##
InstallMethod(Draw,
  "a canvas",
  true,
  [IsCanvas],
  0,
function(canvas)
  local object;
  object := rec();
  object!.mime    := FrancyMIMEType;
  object!.version := InstalledPackageVersion("francy");
  object!.canvas  := Sanitize(canvas);
  return Objectify(
    JupyterRenderableType,
    rec(
      data := rec((FrancyMIMEType) := GapToJsonString(object)),
      metadata := rec((FrancyMIMEType) := rec())
    )
  );
end);

#############################################################################
##
#M  DrawSplash( ) . . . . .
##
InstallMethod(DrawSplash,
  "a canvas",
  true,
  [IsCanvas],
  0,
function(canvas)
    local name, result, page;

    name := Filename(DirectoryTemporary(), Concatenation("francy_", LowercaseString(ReplacedString(canvas!.title, " ", "_")) ,".html"));

    result := Draw(canvas);

    page := Concatenation(
    "<!DOCTYPE html>\n\
    \n\
      \n\
        <meta charset=\"utf-8\" content=\"text/html\" property=\"GAP,francy,d3.v5\"></meta>\n\
        <link rel=\"stylesheet\" type=\"text/css\" href=\"https://cdn.rawgit.com/mcmartins/francy/develop/js/extensions/browser/index.css\"></link>\n\
        \n\
        \n\
        <title>Francy</title>\n\
      \n\
      \n\

\n\
        \n\
      \n\
    ");

    PrintTo(name, page);

    if ARCH_IS_MAC_OS_X() or ARCH_IS_UNIX() then
        Exec("open ",name);
    elif ARCH_IS_WINDOWS() then
        Exec("start ",name);
    fi;

    return page;
end);

max_stars_repo_name: LaGuer/francy

max_issues_repo_name: LaGuer/francy

max_forks_repo_name: LaGuer/francy