
/* 
 * Magma code associated with "Reduction for characters of finite algebra groups," 
 * appearing in J. Algebra 325 (2011), no. 1, 321--351.
 *
 * Author: Anton Evseev
 * Contact: See http://web.mat.bham.ac.uk/A.Evseev/
 *
 */



// Version of 25 July 2010

// This code works in version 2.16-7 of magma, but may not work in earlier versions.

//The code below realises the algorithm described in Secions 3 and 4 of the paper "Reductions for characters of finite algebra groups". Probably, the most useful functions are MTriangularDegrees, DegreesFromZeros and PatternDegrees, which are located at the bottom of the file and are referred to as "end-user functions" in these comments. 

//These comments are not complete; if you require further information, please contact me.

// The pattern algebras encoded by end-user functions are given by the following posets:
// 1. For MTriangularDegrees(n), the poset is the total ordering on n elements.
// 2. For DegreesFromZeros(index,zero), the poset is (index,R) where R is the order obtained from the < order on index by removing all entries of the list zero (it is assumed that index consists of integers).
// 3. For PatternDegrees(ind,C) ,the poset is (ind,R) where R is the transitive closure of C.

// The comments below refer to the notions of "algebraic data" and "categorisation" defined in the paper; in the code these are represented in a different but equivalent way to the manner they are described in the paper.

// Algebraic data A are given by a triple <basis,prod,par> where
// (i) the basis [x_1,...x_N] (the x_i are entities of the same type) such that x_i x_j is a linear combination of elements x_k with k greater than i,j. For technical reasons, each x_i must be of the form <a,b> where a and b are integers (cf. the function outerconjugate_adhoc). 
// (ii) the sequence prod consisting of tuples <i,j,k,c>, where x_k appears in x_i x_j with a non-zero coefficient c. c is given as a sequence of numbered parameters (i.e. elements of par`names) to be multiplied together (so prod encodes the map R of the paper).
// (iii) the list par specifying parameters and relations between them. The field `names' identifies the parameters; `eqns' gives linear relations between them (the first corresponds to 1, the others to parameters in the order of `names'); `nz' lists the parameters that must be non-zero; `mons' is a list of items of the form <i, [j_1,...,j_s]>: the i-th parameter must be the product of parameters j_1,...,j_s. Note that any finite set of polynomial equations and inequalities can be specified in this manner.

// Throughout, hist tracks (in a certain way) the recursive steps that lead to a particular branch of the algorithm and can be ignored for many purposes. 

// Each of the end-user functions finds a categorisation O that is given by a triple res, exc, ec where
// (i) res is a polynomial in two variables x and y; it expresses the characters (or, more precisely, AC-pairs) that the program has been able to "count", i.e. the first entry of O_0 in the notation of the paper. 
// The coefficient in y^m in res is the polynomial F(x), where F(q) is the number of irreducible characters of degree q^m that are being counted (in Section 3 of the paper q and t take the place of x and y)
// (ii) exc is a list [E_1,...E_r] where E_i corresponds to O_i of the paper. As described in the paper, O_i=(F_i,k_i,l_i,m_i), and we have
// E_i=[* A_i, <k_i,l_i,m_i>, z_i, hist_i *] where F_i=Irr(A_i,z_i) and hist_i can largely be ignored (see the previous paragraph).
// (iii) ec consists of entries of the form [* <u_i,v_i,e_i>, par, hist *], each of which corresponds to an entry U_i of O_0 
// (so if O_0 = (f,U_1,...,U_n) then ec has n entries). If U_i=(Q_i,E_i,u_i,v_i,e_i) then par describes (Q_i,E_i) (as explained above), and u_i,v_i,e_i are as given.

// However, the output of "end-user" functions contains slightly more information than a categorisation: the functions also deal with the cases entries E_i of exc such that the algebra encoded by A_i is isomorphic to xF_q[x]/(x^3) (cf. Theorem 2.8 [2.7 in earlier versions] of the paper); let us call such entries E_i "special". The output is of the form res_new, exc_new, ec, ec_spec, res_spec, exc where
// (i) ec and exc are as above;
// (ii) for each special entry E_i the program attempts to compute the total number of characters "counted" by the E_i; this is done by trying to count the number of solutions of equations and inequalities encoded in A_i and, if successful, multiply the outcome by x(x-1) and by (x-1)^{k_i}x^{l_i}y^{m_i}. The polynomial res_spec in variables x and y encodes the numbers of degrees of such "special" characters, for which the count was successful. The data ec_spec is a list of the cases where the count has failed; for each such E_i, the corresponding entry of ec_spec contains the par entry of A_i and <k_i+1,l_i+1,m_i> instead of <k_i,l_i,m_i> (to account for the fact that 1+xF_q[x]/(x^3) has x(x-1) characters whose kernel does not contain <z>, where z is the basis vector given by x^2).
// (iii) res_new = res+res_spec;
// (iv) exc_new is obtained from exc by removing all special entries E_i;

// Thus, if one wishes to count all of the original pattern algebra and sort them by degrees, one should find a way to count all characters encoded by exc_new, by ec and by ec_spec and add these to res_new.



MV<x,y>:=PolynomialRing(Integers(),2);

algtype:=recformat<basis, prod, par>;
partype:=recformat<names,eqns,nz,mons>;
partype_alt:=recformat<names,eqns,nz,rels>;
pattalg:=recformat<index,basis>;

// Splits the algebra A into cases. In each case it is determined which structure constants are zero and which are not;
// This function needs tidying up

split_into_cases:=function(A,hist)
	for T in A`prod do
		for coeff in T[4] do if not coeff in A`par`nz then
			B:=rec<algtype | basis:=A`basis, prod:=A`prod, par:=A`par>;
	 		Append(~B`par`nz,coeff);
			Bcons:=true;
		// checking B for consistency
			j:=Position(B`par`names,coeff);
			for eqn in B`par`eqns do
				for i:=1 to #B`par`names+1 do
					if not i eq j+1 and not eqn[i] eq 0 then continue eqn; end if; 
				end for;
				if not eqn[j+1] eq 0 then Bcons:=false; break eqn; end if;
			end for;
			zeros:=[coeff];
			C:=rec<algtype | basis:=A`basis, prod:=[], par:=A`par>;
			repeat
				change:=false;
				MonsToRemove:=[];
				for zco in zeros do for monno:=1 to #C`par`mons do if zco in C`par`mons[monno][2] then 
					if not C`par`mons[monno][1] in zeros then 
						Append(~zeros,C`par`mons[monno][1]); 
						change:=true; 
						Append(~MonsToRemove,monno);
					end if;
				end if; end for; end for;
				for i:=#MonsToRemove to 1 by -1 do Remove(~C`par`mons,MonsToRemove[i]); end for;
			until not change;
			Ccons:=true;
			for zco in zeros do
				pos:=Position(C`par`names,zco);
				rem:=true;
				for mon in C`par`mons do 
					if mon[1] eq zco then
						// adding the equation saying that zco=0
						eqn:=[0: i in [1..#C`par`names+1]];
						eqn[pos+1]:=1;
						Append(~C`par`eqns,eqn);
						rem:=false;
						break mon;
					end if;
				end for;
				if rem then    //we can remove the variable zco
					if zco in C`par`nz then Ccons:=false; end if;
					for i:=1 to #C`par`eqns do Remove(~C`par`eqns[i],pos+1); end for; 
					Remove(~C`par`names,pos);
				end if;
			end for;
		// checking C for consistency
			if Ccons then
				for eqn in C`par`eqns do
					nze:=0;
					for i:=1 to #C`par`names+1 do
						if not eqn[i] eq 0 then nze:=i; break i; end if; 
					end for;
					if nze gt 0 and (nze eq 1 or C`par`names[nze-1] in C`par`nz) then 
						for i:=nze+1 to #C`par`names+1 do 
							if not eqn[i] eq 0 then continue eqn; end if;
						end for;
 						Ccons:=false;	
						break eqn;
					end if;
				end for;
			end if;
		//
			for S in A`prod do
				 for c in S[4] do if c in zeros then continue S; end if; end for; 
				 Append(~C`prod,S); 
			end for;
			if Ccons then spC:=$$(C,hist); else spC:=[]; end if;
			if Bcons then spB:=$$(B,hist); else spB:=[]; end if; 
			return spC cat spB; 
		end if; end for;
	end for;
	return [A];
end function;




// The function below is used to try to speed up the algorithm and is not described in the paper (though a similar trick for pattern algebras is).
// If it returns true then the map pair[1] -> pair[1]+a*pair[2] can be extended to an isomorphism of A for any value of the parameter a. (N.B. This is an "if", not an "if and only if" statement.)
//The basis of A must be labeled by symbols <i,j> with i and j integers.

outerconjugate_adhoc:=function(A,pair)
	if pair[1][1] eq pair[2][1] then cp:=<pair[1][2],pair[2][2]>;
		elif pair[1][2] eq pair[2][2] then cp:=<pair[2][1],pair[1][1]>;
		else return false;
	end if;
	if <cp[2],cp[1]> in A`basis or <cp[1],cp[2]> in A`basis then return false; end if;
	Tok:=[];
	for i:=1 to #A`prod do if not i in Tok then
		T:=A`prod[i]; 
		if T[3][2] eq cp[1] and <T[3][1],cp[2]> in A`basis then 
			if not T[1][2] eq T[2][1] or not T[1][1] eq T[3][1] or not T[2][2] eq cp[1] then return false; end if;
			j:=Position(A`prod,<T[1],<T[2][1],cp[2]>,<T[3][1],cp[1]>,T[4]>);
			if j eq 0 or j in Tok then return false;
			else 
				Tok cat:=[i,j]; 
				continue i;
			end if; 
		end if;
		if T[3][1] eq cp[2] and <cp[1],T[3][2]> in A`basis then 
			if not T[1][2] eq T[2][1] or not T[2][2] eq T[3][2] or not T[1][1] eq cp[2] then return false; end if;
			j:=Position(A`prod,<<cp[1],T[1][2]>,T[2],<cp[1],T[3][2]>,T[4]>);
			if j eq 0 or j in Tok then return false;
			else 
				Tok cat:=[i,j]; 
				continue i;
			end if; 
		end if;
		if T[2][2] eq cp[2] and <T[2][1],cp[1]> in A`basis then
			if not T[1][1] eq T[3][1] or not T[1][2] eq T[2][1] or not T[3][2] eq cp[2] then return false; end if;
			j:=Position(A`prod,<T[1],<T[2][1],cp[1]>,<T[3][1],cp[1]>,T[4]>);
			if j eq 0 or j in Tok then return false;
			else 
				Tok cat:=[i,j]; 
				continue i;
			end if;
		end if;
		if T[1][1] eq cp[1] and <cp[2],T[1][2]> in A`basis then
			if not T[1][2] eq T[2][1] or not T[2][2] eq T[3][2] or not T[3][1] eq cp[1] then return false; end if;
			j:=Position(A`prod,<<cp[2],T[1][2]>, T[2], <cp[2],T[3][2]>, T[4]>); 
			if j eq 0 or j in Tok then return false;
			else 
				Tok cat:=[i,j]; 
				continue i;
			end if;
		end if;
		if T[2][1] eq cp[1] and <cp[2],T[2][2]> in A`basis then
			if not T[1][1] eq T[3][1] or not T[1][2] eq cp[1] or not T[2][2] eq T[3][2] then return false; end if;
			j:=Position(A`prod,<<T[1][1],cp[2]>,<cp[2],T[2][2]>,T[3],T[4]>); 
			if j eq 0 or j in Tok then return false;
			else 
				Tok cat:=[i,j]; 
				continue i;
			end if;
		end if;
		if T[1][2] eq cp[2] and <T[1][1],cp[1]> in A`basis then
			if not T[1][1] eq T[3][1] or not T[2][1] eq cp[2] or not T[2][2] eq T[3][2] then return false; end if;
			j:=Position(A`prod,<<T[1][1],cp[1]>,<cp[1],T[2][2]>, T[3],T[4]>);
			if j eq 0 or j in Tok then return false;
			else 
				Tok cat:=[i,j]; 
				continue i;
			end if;
		end if;
	end if; end for;
	return true;
end function;

// All the functions from here until proper_count are part of a routine that attempts to count the number of solutions of a system of polynomial equations and inequations (see the first paragraph of the description of the function General in the paper).

// var must be in P`nz
pDivide:=function(P,var,todiv)
	//print "Divide", P, var, todiv;
	if not var in P`nz then return false, _; end if;
	posns:=[Position(P`names,todiv[i]): i in [1..#todiv]];
	posvar:=Position(P`names,var);
	diveqns:=[];
	repeat
		changed:=false;
		for k:=1 to Nrows(P`eqns) do
			for i in posns do
				if not P`eqns[k, i+1] eq 0 then
					Include(~diveqns,k);
					if not P`eqns[k, 1] eq 0 then return false,_; end if; 
					for j:=1 to #P`names do
						if not P`eqns[k, j+1] eq 0 and not j in posns and not j eq posvar then 
							Append(~posns,j); 
							changed:=true; 
						end if;
					end for; 
				end if;
			end for;
		end for;
		for rel in P`rels do
			cnt:=0;
			varcnt:=0;
			for multiple in rel[1] do
				posmlt:=Position(P`names,multiple);
				if posmlt in posns then cnt-:=1; 
					elif multiple eq var then varcnt-:=1;
				end if;
			end for;
			for multiple in rel[2] do
				posmlt:=Position(P`names,multiple);
				if posmlt in posns then cnt+:=1; 
					elif multiple eq var then varcnt+:=1;
				end if;
			end for;
			if cnt lt 0 and cnt+varcnt lt 0 then s:=2; 
				elif cnt gt 0 and cnt+varcnt gt 0 then s:=1;
				else s:=0;
			end if;
			if not s eq 0 then 
				for multiple in rel[s] do if not Position(P`names,multiple) in posns and not multiple eq var then
					newtodiv:=Append([P`names[i]: i in posns], multiple);
					succ, Q:=$$(P,var,newtodiv);
					if succ then 
						return true, Q; 
					end if;
				end if; end for;
				return false, _; 
			end if;
		end for;
		// if posvar in posns then return false,_; end if;
	until not changed;
	todiv:=[P`names[i]: i in posns];
	simpl:=false;
	for i:=1 to #P`rels do
		cnt:=0;
		for s:=1 to 2 do
			for multiple in P`rels[i][s] do if multiple in todiv then cnt+:=2*s-3; end if; end for;
		end for;
		if cnt lt 0 then 
			for tmp:=1 to -cnt do Exclude(~P`rels[i][2],var); end for;
			simpl:=true;
		elif cnt gt 0 then
			for tmp:=1 to cnt do Exclude(~P`rels[i][1],var); end for;
			simpl:=true;
		end if;
	end for;
	for i in diveqns do
		if not P`eqns[i, posvar+1] eq 0 then
			P`eqns[i, 1]+:=P`eqns[i, posvar+1];
			P`eqns[i, posvar+1]:=0;
			simpl:=true;
		end if;
	end for;
	if simpl then return true, P; else return false,_; end if;
end function;


// used by proper_count_main

pforced_zeros:=function(P,zeros)
	//print "forced_zeros", P, zeros;
	for var in zeros do
		if var in P`nz then return []; end if;
	end for;
	for i:=1 to #P`rels do 
		for s:=1 to 2 do for multiple in P`rels[i][s] do
			if multiple in zeros then
				for a in P`rels[i][3-s] do if a in zeros then 
					Remove(~P`rels,i);
					return $$(P,zeros);
				end if; end for;
				seq:=[];
				for a in P`rels[i][3-s] do if not a in P`nz then 
					Q:=P;
					zerosQ:=Append(zeros,a);
					Remove(~Q`rels,i);
					seq cat:=$$(Q,zerosQ);
					Append(~P`nz, a);
				end if; end for;
				return seq;
			end if;
		end for; end for;
	end for;
	for a in zeros do if a in P`nz then return []; end if; end for;
	nos:=[Position(P`names,zeros[i]): i in [1..#zeros]];
	P`eqns:=Submatrix(P`eqns, [1..Nrows(P`eqns)], [1] cat [i+1: i in [1..#P`names] | not i in nos]);
	P`names:=[P`names[k]: k in [1..#P`names] | not k in nos];
	return [P];
end function;


pcleanup:=function(P)
	//print "cleanup", P;
	succ:=-1; //no change
	P`eqns:=HermiteForm(P`eqns);
	repeat
		//print "-------------------------------", P;
		changed:=false;
		donesmth:=false;
		specchange:=false;
		oldP:=P;
		for i:=1 to #P`rels do
			if #P`rels[i][1] ge 2 or #P`rels[i][2] ge 2 then continue i; 
			elif #P`rels[i][1] eq 1 and #P`rels[i][2] eq 1 and not P`rels[i][1][1] eq P`rels[i][2][1] then
					var:=P`rels[i][1][1];
					repvar:=P`rels[i][2][1];
					pos:=Position(P`names,var);
					reppos:=Position(P`names,repvar);
					sign:=P`rels[i][3];
			elif #P`rels[i][1] eq 1 and #P`rels[i][2] eq 0 then
				var:=P`rels[i][1][1];
				pos:=Position(P`names,var);
				reppos:=0;
				sign:=P`rels[i][3];
			elif #P`rels[i][2] eq 1 and #P`rels[i][1] eq 0 then
				var:=P`rels[i][2][1];
				pos:=Position(P`names,var);
				reppos:=0;
				sign:=P`rels[i][3];
			else 
				if P`rels[i][3] eq -1 then continue i; end if;
				donesmth:=true;
			end if;
			Remove(~P`rels,i);	
			if not donesmth then changed:=true; end if;
			break i;
		end for;
		if not changed then 
			for i:=1 to Nrows(P`eqns) do
				cnt:=0; 
				for k:=1 to #P`names+1 do if not P`eqns[i, k] eq 0 then
					if not P`eqns[i, k] eq 1 and not P`eqns[i, k] eq -1 then continue i; end if;
					cnt+:=1;
					if cnt eq 1 then
						reppos:=k-1;
						if reppos gt 0 then repvar:=P`names[reppos]; end if;
						sign:=-P`eqns[i][k];
					elif cnt eq 2 then
						pos:=k-1;
						var:=P`names[pos];
						sign*:=P`eqns[i][k];
					end if;
				end if; end for;
				if cnt eq 2 then 
					changed:=true;
					RemoveRow(~P`eqns,i);
					break i;
				end if;
			end for;
		end if;
		if not changed then 
			for i:=1 to #P`rels do 
				for a in P`rels[i][1] do 
					if a in P`nz and a in P`rels[i][2] then
						Exclude(~P`rels[i][1],a);
						Exclude(~P`rels[i][2],a);
						donesmth:=true;
					end if;
				end for;
			end for;
		end if;
		if changed then
			succ:=1;
			if var in P`nz and reppos gt 0 then Include(~P`nz,repvar); end if;
			for j:=1 to Nrows(P`eqns) do
				P`eqns[j, reppos+1]+:=sign*P`eqns[j][pos+1];
			end for;
			RemoveColumn(~P`eqns,pos+1);
			if reppos eq 0 then
				for i:=1 to #P`rels do
					for s:=1 to 2 do 
						for l:=1 to #P`rels[i][s] do
							if P`rels[i][s][l] eq var then P`rels[i][3]*:=sign; end if;
						end for;
						P`rels[i][s]:=[a: a in P`rels[i][s] | not a eq var];
					end for;
				end for;
			else
				for i:=1 to #P`rels do
					for s:=1 to 2 do 
						for l:=1 to #P`rels[i][s] do
							if P`rels[i][s][l] eq var then 
								P`rels[i][3]*:=sign; 
								P`rels[i][s][l]:=repvar;
							end if;
						end for;
					end for;
				end for;
			end if;
			Remove(~P`names,pos);
			Exclude(~P`nz,var);
		end if;
		if not changed and not donesmth then
			for i:=1 to #P`rels do
				if #P`rels[i][1] eq 1 then s:=1; 
					elif #P`rels[i][2] eq 1 then s:=2;
					else continue i;
				end if;
				var:=P`rels[i][s][1];
				posvar:=Position(P`names,var);
				for k:=1 to Nrows(P`eqns) do if not P`eqns[k,posvar+1] eq 0 then continue i; end if; end for;
				repvars:=P`rels[i][3-s];
				if var in repvars then continue i; end if;
				sign:=P`rels[i][3];
				changed:=true;
				no:=i;
				break i;
			end for;
			if changed then
				succ:=1;
				if var in P`nz then 
					for a in repvars do Include(~P`nz,a); end for; 
					Exclude(~P`nz,var);
				end if;
				Remove(~P`rels,no);
				Remove(~P`names,posvar);
				RemoveColumn(~P`eqns,posvar+1);
				for i:=1 to #P`rels do for s:=1 to 2 do
					tmpcnt:=#P`rels[i][s];
					toremove:=[];
					for j:=1 to tmpcnt do if P`rels[i][s][j] eq var then
						Append(~toremove,j);
						P`rels[i][s] cat:=repvars;
						P`rels[i][3]*:=sign;
					end if; end for;
					P`rels[i][s]:=[P`rels[i][s][j]: j in [1..#P`rels[i][s]] | not j in toremove];
				end for; end for;
			end if;
		end if;
		if not changed and not donesmth then
			for i:=1 to #P`rels do for j:=1 to #P`rels do if not i eq j then
				for s1:=1 to 2 do for s2:=1 to 2 do 
					if SequenceToMultiset(P`rels[i][s1]) subset SequenceToMultiset(P`rels[j][s2]) and #P`rels[i][3-s1] lt #P`rels[i][s1] then
						// replacing P`rels[i][s1] with P`rels[i][3-s1] in P`rels[j][s2]
						for k:=1 to #P`rels[i][s1] do Exclude(~P`rels[j][s2],P`rels[i][s1][k]); end for;
						P`rels[j][s2] cat:=P`rels[i][3-s1];
						P`rels[j][3]*:=P`rels[i][3];
						changed:=true;
						succ:=1;
						break i;
					end if;
				end for; end for;
			end if; end for; end for;
		end if;
		//if changed or donesmth then print P; end if;
	until not changed and not donesmth and not specchange;
	return succ, P;
end function;

pNegate:=procedure(~comps)
	comps[1]:=-comps[1];
	for i:=2 to #comps do comps[i]:=[* comps[i], -1 *]; end for;
end procedure;

pAddComp:=procedure(~comps,newcomp)
	comps[1]+:=newcomp[1];
	comps cat:=[* newcomp[i]: i in [2..#newcomp] *];
end procedure;



proper_count_main:=function(P)
	//print "proper_count_main", P;
	if Nrows(P`eqns) eq 0 and #P`rels eq 0 then return true, (x-1)^(#P`nz)*x^(#P`names-#P`nz); end if;
	//succ, Q:=pcleanup(P);
	//if succ eq 1 then P:=Q; end if;
	comps:=[* 0 *];  //data to be returned if unsuccessful
	lst:=[i: i in [1..#P`names] | not P`names[i] in P`nz] cat [i: i in [1..#P`names] | P`names[i] in P`nz];
// finding a variable that appears only once
	for i:=1 to #P`names do
		no:=0;
		for j:=1 to Nrows(P`eqns) do if not P`eqns[j,i+1] eq 0 then 
			if not P`eqns[j,i+1] eq 1 and not P`eqns[j,i+1] eq -1 then continue i; end if; // ! see original
			if no gt 0 then continue i; end if;
			no:=1;
			type_app:=1;
			eqno:=j;
		end if; end for;
		for j:=1 to #P`rels do
			pos:=Position(P`rels[j][1],P`names[i]);
			if pos gt 0 then
				if no gt 0 then continue i; end if;
				if Position(P`rels[j][1],P`names[i],pos+1) gt 0 then continue i; end if;  
				no:=1;
				type_app:=2;
				part:=1;
				relno:=j;
			end if;
			pos:=Position(P`rels[j][2],P`names[i]);
			if pos gt 0 then
				if no gt 0 then continue i; end if;
				if Position(P`rels[j][2],P`names[i],pos+1) gt 0 then continue i; end if;
				no:=1;
				type_app:=2;
				part:=2;
				relno:=j;
			end if;
		end for;
		if no eq 1 then 
			case type_app:
			when 1:
				if P`names[i] in P`nz then
					cnt:=0;
					other:=0;
					for k:=1 to #P`names+1 do if not P`eqns[eqno,k] eq 0 then 
						cnt+:=1; 
						if not k eq i+1 then other:=k; end if;	
					end if; end for;
					Exclude(~P`nz,P`names[i]);
					if cnt le 2 then
						if other eq 0 then return true, 0; end if;
						RemoveRow(~P`eqns,eqno);
						RemoveColumn(~P`eqns,i+1);
						if other gt 1 then Include(~P`nz,P`names[other-1]); end if;
						Remove(~P`names,i);
						return $$(P);
					end if;
					if cnt gt 2 then 
						//print P, i;
						seq:=pforced_zeros(P,[P`names[i]]);
						sumres:=0; 
						totalsucc:=true;
						for R in seq do
							succ,res:=$$(R);
							if succ then sumres-:=res;
							else
								pNegate(~res);
								pAddComp(~comps,res);
								totalsucc:=false;
							end if;
						end for;
						RemoveRow(~P`eqns,eqno);
						RemoveColumn(~P`eqns,i+1);
						Remove(~P`names,i);
						succ,nzres:=$$(P);
						if succ and totalsucc then return true, nzres+sumres;  
						elif not succ then
							comps[1]+:=sumres;
							pAddComp(~comps,nzres);
							return false, comps;
						else 
							comps[1]+:=sumres+nzres;
							return false,comps;
						end if;
					end if;
				else
					RemoveRow(~P`eqns,eqno);
					RemoveColumn(~P`eqns,i+1);
					Remove(~P`names,i);
					return $$(P);
				end if;
			when 2:
				var:=P`names[i];
				R:=P;
				sumres:=0;
				totalsucc:=true;
				for multiple in R`rels[relno][part] do if not multiple eq var and not multiple in R`nz then
					// the case when multiple = 0
					seqR:=pforced_zeros(R,[multiple]);
					for Q in seqR do
						succ, res:=$$(Q);
						if succ then sumres+:=res; else 
							pAddComp(~comps,res); 
							totalsucc:=false; // !
						end if;
					end for;
						Append(~R`nz,multiple);
					end if; 
				end for;
				Remove(~R`names,i);
				RemoveColumn(~R`eqns,i+1);
				nzpos:=Position(R`nz,var);
				if nzpos gt 0 then 
					Remove(~R`nz,nzpos); // ! redundant "Exclude" was here
					for a in R`rels[relno][3-part] do Include(~R`nz,a); end for;
				end if;
				Remove(~R`rels,relno);
				succ, nzres:=$$(R);
				if succ and totalsucc then return true, nzres+sumres; 
				elif not succ then 
					comps[1]+:=sumres;
					pAddComp(~comps,nzres); 
					return false, comps;
				else 	comps[1]+:=sumres+nzres;
					return false,comps; 
				end if;
			end case;
		end if;
	end for;	
// trying to get rid of relations by a simplification 
	if #P`rels gt 0 then 
		for rel in P`rels do for s in [1,2] do 
			for var in rel[s] do if var in P`nz then
				succ,Q:=pDivide(P,var,rel[3-s]);
				// if succ then print "________________", P, var, mon[1], Q; end if;
				if succ then
					return $$(Q); 
				end if;
			end if; end for;
		end for; end for;
		succ, Q:=pcleanup(P);
		if succ eq 1 then 
			return $$(Q); 
		else return false, [* 0, P *]; end if;
	end if; 
// assume that there aren't any monomials
	PS:=[SetToSequence(Su): Su in Subsets({i: i in [1..#P`names]| P`names[i] in P`nz})];
	res:=0;
	for T in PS do
		// new code here 16.07.2009
		B:=Submatrix(P`eqns, [1..Nrows(P`eqns)],[i: i in [1..#P`names+1] | not i-1 in T]);
		//
		eldiv:=ElementaryDivisors(B);
		for d in eldiv do if not d eq 1 then print "*******************", P, B; return false, [* 0, P *]; end if; end for;
		C:=RemoveColumn(B,1);
		eldivC:=ElementaryDivisors(C);
		for d in eldivC do if not d eq 1 then print "******************", P; return false, [* 0, P *]; end if; end for;
		if #eldivC eq #eldiv then res+:=(-1)^(#T)*x^(#P`names-#T-#eldivC); end if;
	end for;
	return true,res;
end function;


proper_count:=function(P)
	Qeqns:=Matrix(Integers(), #P`eqns, #P`names+1,P`eqns);
	Qrels:=[ <[mon[1]], mon[2], 1> : mon in P`mons ];
	Q:=rec<partype_alt | names:=P`names, eqns:=Qeqns, nz:=P`nz, rels:=Qrels>;
	return proper_count_main(Q);
end function;



forward degspec;

// degnonspec is the function General of the paper.
// It outputs a categorisation of the character degrees of an algebra A and multiplies the output by the polynomial f (working recursively).
// f is be given as a triple <k,l,m> corresponding to the polynomial (x-1)^k x^l y^m; when we say that the output is multiply by f, we mean that all the character counts are multiplied by (x-1)^k*x^l and all the character degrees are multiplied by y^m. 

// Recall: all elements of A`basis must be of the form <i,j> where i and j are integers.

degnonspec:=function(A,f: hist:=[* [], [], [], [] *])
	if #A`prod eq 0 then 
		succ,res:=proper_count(A`par);
		//if #A`par`eqns gt 0 or #A`par`mons gt 0 then print "______________________________", A`par, res; end if;
		if succ then
			return (x-1)^f[1]*x^f[2]*y^f[3]*res*x^#A`basis, [], []; 
		else return 0, [], [ [* <f[1],f[2]+#A`basis,f[3]>, A`par, hist *] ];
		end if; 
	end if;
// finding a good N
	Z:=A`basis;
	for T in A`prod do Exclude(~Z,T[1]); Exclude(~Z,T[2]); end for;
	N:=Z[#Z];
	CanBeRemoved:=[ []: i in [1..#Z]];
	for i:=1 to #Z do
		for K in Z do if not K eq Z[i] and outerconjugate_adhoc(A,<K,Z[i]>) then 
			Append(~CanBeRemoved[i],K); 
	//test
	//		if not outerconjugate(A,<K,Z[i]>) then print "***********************", A, K, Z[i]; end if;
		end if; end for;
	end for;
	min:=-1;
	maxrem:=0;
	for possNno:=#Z to 1 by -1 do
		leftmult:=[];
		count:=0;
		for T in A`prod do if T[3] eq Z[possNno] then Include(~leftmult,T[1]); end if; end for;
		for T in A`prod do if T[2] in leftmult then count+:=1; end if; end for;
		// if count eq 0 then N:=Z[possNno]; break; end if;
		if min eq -1 or count lt min or (count eq min and #CanBeRemoved[possNno] gt maxrem) then
			min:=count; 
			maxrem:=#CanBeRemoved[possNno]; 
			Nno:=possNno; 
		end if; 
	end for;
	N:=Z[Nno];
	C:=A;
	C`basis:=[A`basis[i]: i in [1..#A`basis] | not A`basis[i] in CanBeRemoved[Nno]];
	C`prod:=[A`prod[i]: i in [1..#A`prod] | not A`prod[i][3] in CanBeRemoved[Nno]];
	fnew:=f; 
	fnew[2]+:=#CanBeRemoved[Nno];
// applying recursion
	// tracking history
	histnew:=hist;
	for K in CanBeRemoved[Nno] do Append(~histnew[3], <N,K,#histnew[1]+#histnew[2]+#histnew[3]+#histnew[4]+1>); end for;
	res1,exc1,exc_count1:=degspec(C,fnew,N: hist:=histnew);
	B:=rec<algtype | basis:=Exclude(A`basis,N), prod:=[], par:=A`par>;
	for T in A`prod do if not T[3] eq N then Append(~B`prod, T); end if; end for;
	histnew:=hist;
	Append(~histnew[4],<N,#hist[1]+#hist[2]+#hist[3]+#hist[4]+1>);
	res2,exc2,exc_count2:=$$(B,f: hist:=histnew);
//test	
//	succ,npar:=count_parameters(A`par);
//	if #exc1 eq 0 and #exc2 eq 0 and #exc_count1 eq 0 and #exc_count2 eq 0 and not MV!Evaluate(res1+res2,[x,x^2]) eq MV !  x^#A`basis*Evaluate((x-1)^f[1]*x^f[2]*y^f[3],[x,x^2])*npar then 
//		print "***********************", A, N;
//	end if;
// end of test	
// MV! Evaluate(res1+res2,[x,x^2]) div (Evaluate(MV! (x-1)^f[1]*x^f[2]*y^f[3],[x,x^2])*npar); 
	return res1+res2, exc1 cat exc2, exc_count1 cat exc_count2;	 
end function;

// used when new parameters are introduced
newindex:=function(List)
	r:=0;
	for i in List do if i gt r then r:=i; end if; end for; 
	r+:=1;
	return r;
end function;

// degspec is the function TypeB of the paper
//finds the degrees of characters \chi of 1+A such that \chi restricted to 1+<x_N> is non-trivial (again, special cases are excepted and returned separately)

// Recall: all elements of A`basis must be of the form <i,j> where i and j are integers.

degspec:=function(A,f,N: hist:=[* [], [], [], [] *])
	mult:=0;
	repeat
		// finding Z, the set of basis vectors that are in Z_0 (A)
		Z:=A`basis;
		for T in A`prod do 
			Exclude(~Z, T[1]);
			Exclude(~Z, T[2]);
		end for;
		// factoring out basic vectors x_K in Z with x_K ->x_K +a*x_N extendible to an automorphism
		change:=false;
		for K in Z do if not K eq N and outerconjugate_adhoc(A,<K,N>) then 
			mult+:=1;
			change:=true;
			// removing K
			Exclude(~A`basis,K);
			A`prod:=[A`prod[i]: i in [1..#A`prod] | not A`prod[i][3] eq K]; 
		// tracking history
			Append(~hist[3],<N,K,#hist[1]+#hist[2]+#hist[3]+#hist[4]+1>);
		end if; end for;
	until not change;
	f[2]+:=mult;
//
	Ms:=[];
	countMs:=[];
	notinMs:=[];
	// finding a two-dimensional ideal <x_M,x_N> with Ax_M=0;
	sep:=true; // sep will record whether <x_1,...,x_{N-1}> is a subalgebra of A.
//
	for T in A`prod do 
		Include(~notinMs,T[2]);
		if not T[3] eq N then Include(~notinMs,T[1]);
		else sep:=false; 
		end if;
	end for;
	for T in A`prod do 
		if T[3] eq N and not T[1] in notinMs then 
			i:=Position(Ms,T[1]); 
			if i eq 0 then Append(~Ms,T[1]); Append(~countMs,1); 
			else countMs[i]+:=1; 
			end if;
		end if;
	end for;
	if #Ms gt 0 then
		temp,i:=Min(countMs);
		M:=Ms[i]; // we won't use the rest of Ms 
		R:=[]; Rcoeffs:=[];
		for T in A`prod do
			if T[3] eq N and T[1] eq M then Append(~R,T[2]); Append(~Rcoeffs,T[4]); end if;
		end for;
// finding the last element in R w.r.t. the given ordering of the basis
		rmpos:=0;
		for k in R do 
			tmppos:=Position(A`basis,k);
			if tmppos gt rmpos then rmpos:=tmppos; end if;
		end for;
		k:=A`basis[rmpos];
// constructing the new algebra B:=C_A (x_M)/<x_M>;
		B:=rec<algtype | basis:=Exclude(A`basis,k), prod:=[], par:=A`par>;
		Exclude(~B`basis,M);
		extras:=[];
		posinR:=Position(R,k);
		ck:=Rcoeffs[posinR];
		Remove(~R,posinR);
		Remove(~Rcoeffs,posinR);
// finding the structure constants (we factor <x_M> out)
		for T in A`prod do if not T[3] eq M and not T[1] eq M and not T[2] eq M and not T[3] eq k then 
			if not T[1] eq k and not T[2] eq k then
				S:=T;
				if T[1] in R then S[4] cat:=ck; end if;
				if T[2] in R then S[4] cat:=ck; end if;
				j:=Position(R,T[3]);
				// dividing by ck
				if not j eq 0 then
					ck_rem:=[];
					for a in ck do 
						i:=Position(S[4],a);
						if not i eq 0 then Remove(~S[4],i); 
							else Append(~ck_rem,a); 	 
						end if;
					end for;
					if #ck_rem gt 0 then
						newnum:=newindex(B`par`names);
						b:=newnum;
						Append(~B`par`names,newnum);
						for k:=1 to #B`par`eqns do Append(~B`par`eqns[k],0); end for;
						newnum+:=1;
						if #S[4] eq 1 then 
							Append(~B`par`mons,<S[4][1], Append(ck_rem,b)>);
						else	// theoretically, S[4] may be empty
							c:=newnum;
							Append(~B`par`names,newnum);
							for k:=1 to #B`par`eqns do Append(~B`par`eqns[k],0); end for;
							newnum+:=1;
							Append(~B`par`mons,<c, S[4]>);
							Append(~B`par`mons,<c, Append(ck_rem,b)>);
						end if;
						S[4]:=[Integers() | b];
					end if;
				end if;
				Append(~B`prod,S);
			elif T[1] eq k then
				if T[2] in R then
					for j:=1 to #R do
						Append(~extras,<R[j],T[2],T[3],T[4] cat Rcoeffs[j] cat ck,-1>);
					end for;
				elif T[2] eq k then
					for i:=1 to #R do for j:=1 to #R do 
						Append(~extras,<R[i],R[j],T[3],T[4] cat Rcoeffs[i] cat Rcoeffs[j],1>);
					end for; end for;
				else  
					for j:=1 to #R do Append(~extras,<R[j],T[2],T[3],T[4] cat Rcoeffs[j],-1>); end for;
				end if;
			elif T[2] eq k then
				if T[1] in R then
					for j:=1 to #R do
						Append(~extras,<T[1],R[j],T[3],T[4] cat Rcoeffs[j] cat ck,-1>);
					end for;
				elif not T[1] eq k then 
					for j:=1 to #R do Append(~extras,<T[1],R[j],T[3],T[4] cat Rcoeffs[j],-1>); end for;
				end if;
			end if;
		end if; end for;
		// dealing with "extras", which encode structure constants to be collated
		extrasexist:=false;
		while #extras gt 0 do
			extrasexist:=true;
			Q:=extras[1];
			already:=false;
			tmp:=[<Q[4],Q[5]>];
			tmpind:=[1];
			for i:=1 to #B`prod do
				if B`prod[i][1] eq Q[1] and B`prod[i][2] eq Q[2] and B`prod[i][3] eq Q[3] then
					already:=true;
					Append(~tmp,<B`prod[i][4],1>);
					ind_inBprod:=i;
					break i;
				end if;	
			end for; 
			for i:=2 to #extras do
				S:=extras[i];
				if S[1] eq Q[1] and S[2] eq Q[2] and S[3] eq Q[3] then
					Append(~tmp,<S[4],S[5]>);
					Append(~tmpind,i);
				end if;	
			end for;
			extras:=[extras[i] : i in [1..#extras] | not i in tmpind];
// collating entries in tmp where possible
			repeat 
				changed:=false; 
				for i:=2 to #tmp do for j:=1 to i-1 do
					if tmp[i][1] eq tmp[j][1] then 
						tmp[j][2]+:=tmp[i][2]; 
						changed:=true;	
						Remove(~tmp,i);
						break i; 
					end if; 
				end for; end for; 
			until not changed;
			tmp:=[tmp[i]: i in [1..#tmp] | not tmp[i][2] eq 0];
//			
			if #tmp eq 0 then
				if already then Remove(~B`prod,ind_inBprod); end if;
			elif #tmp eq 1 and tmp[1][2] eq 1 and not already then
				Append(~B`prod,<Q[1],Q[2],Q[3],Q[4]>); 
			elif #tmp ge 1 then
				newnum:=newindex(B`par`names);
				c:=newnum;
				Append(~B`par`names,newnum);
				for k:=1 to #B`par`eqns do Append(~B`par`eqns[k],0); end for;
				newnum+:=1;
				v:=[0: i in [1..#B`par`names+1]];
				v[#B`par`names+1]:=1;
				for summand in tmp do
					if #summand[1] eq 0 then v[1]-:=summand[2];
					elif #summand[1] eq 1 then v[Position(B`par`names,summand[1][1])+1]-:=summand[2];
					else
						Append(~B`par`names,newnum);
						for k:=1 to #B`par`eqns do Append(~B`par`eqns[k],0); end for;
						Append(~B`par`mons,<newnum,summand[1]>);
						newnum+:=1;
						Append(~v,-summand[2]);
					end if;
				end for;
				Append(~B`par`eqns,v);
				if already then B`prod[ind_inBprod][4]:=[c]; 
				else Append(~B`prod,<Q[1],Q[2],Q[3],[c]>);
				end if;
			end if;
		end while;
// should check whether there could possibly be more than one case first
		cases:=split_into_cases(B, hist);
//		if not succ then return 0, [ [* A,f,N,hist *]], []; end if; // not needed any more
		res:=0;
		exc:=[];
		exc_count:=[];
		B1:=B;
		Append(~hist[1],<N,M, #hist[1]+#hist[2]+#hist[3]+#hist[4]+1>);
		for B in cases do
			fnew:=f;
			fnew[3]+:=1;
			res1,exc1,exc_count1:=$$(B,fnew,N: hist:=hist);
			res+:=res1;
			exc cat:=exc1;
			exc_count cat:=exc_count1;
		end for; 
		return res,exc,exc_count;
	end if;
	if sep then
		B:=rec<algtype | basis:=Exclude(A`basis,N), prod:=A`prod, par:=A`par>;
		fnew:=f;
		fnew[1]+:=1;
		return degnonspec(B,fnew: hist:=hist);
	end if;	
//
// finding a potential M and the central vectors that need to be "factored out" (these are recorded in "quot").
//
	cand:=A`basis;
	for T in A`prod do 		
		Exclude(~cand,T[2]);
		if not T[3] in Z then 
			Exclude(~cand,T[1]);
		end if; 
	end for;
	cand:=[cand[i]: i in [1..#cand] | not cand[i] in Z];
	candspec:=[ [] : i in [1..#cand]];
	for T in A`prod do 
		i:=Position(cand,T[1]); 
		if i gt 0 then Include(~candspec[i],T[3]); end if;
	end for;
// we make sure x_N is involved in x_M*A, if possible	
	Good:=[i: i in [1..#cand] | N in candspec[i] ];	
	Mpos:=0; minival:=0;
	if #Good gt 0 then 
		for i:=1 to #candspec do 
			if i in Good and (Mpos eq 0 or minival gt #candspec[i]) then Mpos:=i; minival:=#candspec[i]; end if;
		end for;
	else
		for i:=1 to #candspec do 
			if Mpos eq 0 or minival gt #candspec[i] then Mpos:=i; minival:=#candspec[i]; end if;
		end for;
	end if;
	if Mpos eq 0 then return 0, [ [* A,f,N,hist *]], []; end if;
	quot:=Exclude(candspec[Mpos],N);
	if #quot eq 0 then print "Error!", A; return 0; end if; // this shouldn't happen
//
// constructing the relations of the new algebra
	newnum:=newindex(A`par`names);
	corr:=AssociativeArray();
	B:=rec<algtype | basis:=A`basis, prod:=[], par:=A`par>;
	for S in quot do 
		Append(~B`par`names,newnum);
	// tracking history
		Append(~hist[2], <N, S, newnum, #hist[1]+#hist[2]+#hist[3]+#hist[4]+1>);
		for k:=1 to #B`par`eqns do Append(~B`par`eqns[k],0); end for;
		corr[S]:=newnum;
		newnum+:=1;
		Exclude(~B`basis,S);
	end for;
	// reordering the basis
	pos:=Position(B`basis,N);
	Remove(~B`basis,pos);
	Append(~B`basis,N);
	//
	tmpind:=[];
	tmppar:=[];
	for T in A`prod do
		if T[3] in quot or T[3] eq N then
			if T[3] eq N then coeff:=T[4]; else coeff:=Append(T[4],corr[T[3]]); end if;
			already:=false;
			for i:=1 to #B`prod do if B`prod[i][1] eq T[1] and B`prod[i][2] eq T[2] and B`prod[i][3] eq N then
				already:=true;
				if #coeff eq 0 then c:=0;
				elif #coeff eq 1 then c:=coeff[1];
				else 
					c:=newnum;
					Append(~B`par`names,newnum);
					for k:=1 to #B`par`eqns do Append(~B`par`eqns[k],0); end for;
					Append(~B`par`mons,<newnum,coeff>);
					newnum+:=1;
				end if;
				pos:=Position(tmpind,<T[1],T[2]>);
				if pos gt 0 then Append(~tmppar[pos],c); 
				else 
					Append(~tmpind,<T[1],T[2]>);
					if #B`prod[i][4] eq 0 then Append(~tmppar,[0,c]);
					elif #B`prod[i][4] eq 1 then Append(~tmppar,[B`prod[i][4][1],c]);
					else	
						Append(~B`par`names,newnum);
						for k:=1 to #B`par`eqns do Append(~B`par`eqns[k],0); end for;
						Append(~B`par`mons,<newnum,B`prod[i][4]>);
						Append(~tmppar,[newnum,c]);
						newnum+:=1;
					end if;	
				end if;
				break i;
			end if; end for;
			if not already then 
				Append(~B`prod,<T[1],T[2],N,coeff>);
			end if;
		else Append(~B`prod,T);
		end if;
	end for; 
	for i:=1 to #tmpind do
		Append(~B`par`names,newnum);
		for k:=1 to #B`par`eqns do Append(~B`par`eqns[k],0); end for;
		inBprod:=false;
		for j:=1 to #B`prod do if tmpind[i][1] eq B`prod[j][1] and tmpind[i][2] eq B`prod[j][2] and B`prod[j][3] eq N then
			B`prod[j][4]:=[newnum]; 
			inBprod:=true;
			break j;
		end if; end for;
		if not inBprod then Append(~B`prod,<tmpind[i][1],tmpind[i][2],N,[newnum]>); end if;
		v:=[0: j in [1..#B`par`names+1]];
		v[#B`par`names+1]:=1;
		for j in tmppar[i] do if j eq 0 then v[1]-:=1; else v[Position(B`par`names,j)+1]-:=1; end if; end for;
		Append(~B`par`eqns,v);
		newnum+:=1;
	end for;
	i:=Position(B`basis,N);
	if not i eq #B`basis then
		Remove(~B`basis,i);
		Append(~B`basis,N);
	end if;
	cases:=split_into_cases(B, hist);
	res:=0;
	exc:=[];
	exc_count:=[];
	for B in cases do
		res1,exc1,exc_count1:=$$(B,f,N: hist:=hist);
		res+:=res1;
		exc cat:=exc1;
		exc_count cat:=exc_count1;
	end for; 
	return res,exc,exc_count;
end function;



// Outputs the character degrees of a pattern algebra P using first "fast reduction" and then, if necessary, degnonspec.
// The basis of a pattern algebra must satisfy the following property: if i\le k\le l \le j and <i,j>,<k,l> are both in
// P`basis then <k,l> comes before <i,j> in the ordering of P`basis.



SortPairsTmp:=function(x,y,B)
	x1B:=Position(B,x[1]);
	y1B:=Position(B,y[1]);
	if x1B lt y1B then return 1; 
		elif y1B lt x1B then return -1;
		else return Position(B,x[2])-Position(B,y[2]);
	end if;
end function;

forward EmptyAlgebra;

// pattern_fr is the function PatternAlgebra described in Section 4 of the paper
// P encodes a poset (index,basis) [so basis is a partial order on index] and f is a triple as in functions degnonspec and defspec. As in those functions, the output is multiplied by the polynomial given by f.

//All elements of P`index must be integers.

pattern_fr:=function(P,f)
	if #P`basis eq 0 then return (x-1)^f[1]*x^f[2]*y^f[3],[],[]; end if;
	v:=P`index[1];
	Q:=rec<pattalg | index:=[P`index[i]: i in [2..#P`index]], basis:=[pair: pair in P`basis | not pair[1] eq v]>;
// trying to speed up
	repeat
		MappedTov0:=[u: u in Q`index | <v,u> in P`basis];
		sep:=0;
		succ_sep:=false;
		for i:=1 to #MappedTov0 do
			for u in Q`index do if <MappedTov0[i],u> in Q`basis or <u,MappedTov0[i]> in Q`basis then continue i; end if; end for;
			sep:=i;
			break;
		end for;
		if not sep eq 0 then	
			succ_sep:=true;
			sep:=MappedTov0[sep];
			f[2]+:=1;
			Exclude(~P`index,sep);
			Exclude(~P`basis,<v,sep>);
			Exclude(~Q`index,sep);
		end if;
	until not succ_sep;
//
// finding (part of) the outer normaliser of Q
	Norm:=[];
	for i:=1 to #Q`index do for j:=i+1 to #Q`index do
		u:=Q`index[i];
		w:=Q`index[j];
		if not <u,w> in Q`basis then
			innorm:=true;
			for k:=1 to i-1 do if not <Q`index[k],w> in Q`basis and <Q`index[k],u> in Q`basis then innorm:=false; end if; end for; 
			for l:=j+1 to #Q`index do if not <u,Q`index[l]> in Q`basis and <w,Q`index[l]> in Q`basis then innorm:=false; end if; end for; 
			if innorm then
				Append(~Norm,<u,w>);
			else
				for k:=1 to j-1 do if not <Q`index[k],u> in Q`basis and <Q`index[k],w> in Q`basis then continue j; end if; end for; 
				for l:=i+1 to #Q`index do if not <w,Q`index[l]> in Q`basis and <u,Q`index[l]> in Q`basis then continue j; end if; end for; 
				Append(~Norm,<w,u>); 
			end if;
		end if;
	end for; end for;
//
	for i:=1 to #Q`index do if <v,Q`index[i]> in P`basis then 
		for j:=#Q`index to i+1 by -1 do 
			if <v,Q`index[j]> in P`basis and <Q`index[j],Q`index[i]> in Norm then
				tmp:=Q`index[i];
				Q`index[i]:=Q`index[j];
				Q`index[j]:=tmp;
			end if;
		end for; end if; 
	end for;

	MappedTov:=[u: u in Q`index | <v,u> in P`basis];
	nzcoeff:=[];
	finished:=false;
	NormUsed:=[];
	KilledByNorm:=[];
	res:=0;
	exc:=[];
	exc_count:=[];
	repeat
		if #nzcoeff eq 0 then
			if #MappedTov gt 0 then    //count the case when nzcoeff=[Q`index[1]] as well
				fnew:=f;
				fnew[2]+:=1;
				res1,exc1,exc_count1:=$$(Q,fnew);
			else
				res1,exc1,exc_count1:=$$(Q,f);
			end if;
		elif #nzcoeff eq 1 then
			fnew:=f;
			R:=Q;
			for u in MappedTov do 
				pos:=Position(R`basis,<u,MappedTov[nzcoeff[1]]>);
				if pos gt 0 then 
					Remove(~R`basis,pos);
					fnew[3]+:=1;
				end if;
			end for;
			fnew[1]+:=1;
			fnew[2]+:=#KilledByNorm;
			res1,exc1,exc_count1:=$$(R,fnew);
		elif #nzcoeff eq 2 then
			v:=MappedTov[nzcoeff[1]];
			u:=MappedTov[nzcoeff[2]];
			A:=EmptyAlgebra();
		//constructing the new basis making sure the ordering satisfies the required condition 
			VectorsToRemove:=[i: i in [1..#Q`basis] | (Q`basis[i][2] eq v and not <Q`basis[i][1],u> in Q`basis and Q`basis[i][1] in MappedTov) or (Q`basis[i][2] eq u and not <Q`basis[i][1],v> in Q`basis and Q`basis[i][1] in MappedTov)];
			A`basis:=[Q`basis[i]: i in [1..#Q`basis] | not i in VectorsToRemove];
			fnew:=f;
			fnew[1]+:=2;
			fnew[2]+:=#KilledByNorm;
			fnew[3]+:=#VectorsToRemove;
			for w in MappedTov do
				pos1:=Position(A`basis,<w,v>);
				if pos1 gt 0 then
					pos2:=Position(A`basis,<w,u>);
					if pos1 lt pos2 then
						Remove(~A`basis,pos2);
					else
						A`basis[pos2]:=<w,v>;
						Remove(~A`basis,pos1);
					end if;
					fnew[3]+:=1;
				end if;
			end for; 
			for kl in A`basis do for lm in A`basis do
				if kl[2] eq lm[1] then
					Append(~A`prod,<kl,lm,<kl[1],lm[2]>,[Integers()| ]>);
				end if;
				if lm[1] eq u and kl[2] eq v and kl[1] in MappedTov then 
					Append(~A`prod,<kl,lm,<kl[1],lm[2]>,[Integers()| ]>); 
				end if;
				if lm[2] eq v and kl[2] eq lm[1] and lm[1] in MappedTov and not kl[1] in MappedTov and <lm[1],u> in Q`basis then
					Append(~A`prod,<kl,lm,<kl[1],u>,[Integers() | ]>);
				end if; 
			end for; end for;
			res1,exc1,exc_count1:=degnonspec(A,fnew: hist:=[* [],[],[],[], A *]);
		else 
			A:=EmptyAlgebra();
			A`basis:=P`basis;	
			for kl in A`basis do for lm in A`basis do
				if kl[2] eq lm[1] then
					Append(~A`prod,<kl,lm,<kl[1],lm[2]>,[Integers()| ]>);
				end if;
			end for; end for;
			return degnonspec(A,f: hist:=[* [],[],[],[],A *]);
		end if;
		res+:=res1;
		exc cat:=exc1;
		exc_count cat:=exc_count1;
		if #nzcoeff eq 0 then
			if #MappedTov le 1 then break; else Append(~nzcoeff,#MappedTov); end if;
			for j:=1 to #MappedTov-1 do 
				if <MappedTov[j],MappedTov[#MappedTov]> in Norm then Append(~KilledByNorm,MappedTov[j]); end if;
			end for; 
		else	
			top:=#MappedTov;
			repeat
				succ:=false;
				if #nzcoeff gt 0 and nzcoeff[#nzcoeff]-1 lt top then top:=nzcoeff[#nzcoeff]-1; end if;
				for i:=top to 1 by -1 do if not MappedTov[i] in KilledByNorm then
					for k in nzcoeff do
						if <MappedTov[i],MappedTov[k]> in Q`basis then continue i; end if;
					end for;
					succ:=true;
					Append(~nzcoeff,i);
					for j:=1 to i-1 do 
						for k:=1 to #nzcoeff-1 do if <MappedTov[j],MappedTov[nzcoeff[k]]> in Q`basis then continue j; end if; end for;
						if <MappedTov[j],MappedTov[i]> in Norm then
							Include(~KilledByNorm,MappedTov[j]);
						end if;
					end for;
					KilledByNorm:=[vect: vect in KilledByNorm | not <vect, MappedTov[i]> in Q`basis];
					break i; 
				end if; end for;
				if not succ then
					top:=nzcoeff[#nzcoeff]-1;
					Prune(~nzcoeff);
					if #nzcoeff eq 0 and top le 1 then finished:=true; break; end if;
					KilledByNorm:=[];
					for j:=1 to #MappedTov-1 do
						killed:=false; 
						for i in nzcoeff do if i gt j then 
							if <MappedTov[j],MappedTov[i]> in Q`basis then continue j; end if;
							if <MappedTov[j],MappedTov[i]> in Norm then killed:=true; end if;
						end if; end for;
						if killed then Append(~KilledByNorm,MappedTov[j]); end if;
					end for;
				end if;
			until succ;
		end if;
	until finished;
// testing
//	if not Evaluate(res,[x,x^2]) eq Evaluate(MV!(x-1)^f[1]*x^f[2]*y^f[3],[x,x^2])*x^#P`basis then print "***********************************", P, Evaluate(res,[x,x^2]) div Evaluate(MV!(x-1)^f[1]*x^f[2]*y^f[3],[x,x^2]); end if;
	return res,exc,exc_count;
end function;

// checks whether the algebras encoded by A are isomorphic to the 2-dimensional commutative algebra xF_q[x]/(x^3)

IsType1:=function(A,N)
	if not #A`basis eq 2 or not A`basis[2] eq N then return false; end if;
	M:=A`basis[1];
	if not #A`prod eq 1 or not A`prod[1][1] eq M or not A`prod[1][2] eq M or not A`prod[1][3] eq N then return false; end if;
	for multiple in A`prod[1][4] do if not multiple in A`par`nz then return false; end if; end for;
	return true;
end function;

// used by AddExceptions

TreatExceptions:=function(exc)
	excres:=0;
	excexc:=[];
	excec:=[];
	for C in exc do
		A:=C[1];
		f:=C[2];
		N:=C[3];
		if IsType1(A,N) then
			succ,cnt:=proper_count(A`par); 
			if succ then excres+:=(x-1)^f[1]*x^f[2]*y^f[3]*cnt*x*(x-1); 
			else
				Append(~excec, [* <f[1]+1,f[2]+1,f[3]>, A`par, C[4] *]);	
			end if;
		else
			Append(~excexc,C);
		end if;
	end for;
	return excres,excexc,excec;
end function;

// deals with core algebras of the form xF_q[x]/(x^3).

AddExceptions:=function(res,exc,ec)
	excres,excexc,excec:=TreatExceptions(exc);
	return res+excres, excexc, ec, excec, excres, exc;
end function;


// creates the algebra of nilpotent triangular n-by-n matrices.
TriangularAlgebra:=function(n)
	A:=rec<algtype | basis:=[], prod:=[], par:=rec<partype| names:=[],eqns:=[],nz:=[],mons:=[]> >;	
	for i:=n to 1 by -1 do for k:=i+1 to n do 
		Append(~A`basis,<i,k>);
	end for; end for;
	for i:=1 to n do for j:=i+1 to n do for k:=j+1 to n do
		Append(~A`prod,<<i,j>,<j,k>,<i,k>,[Integers() | ]>);
	end for; end for; end for;
	return A;
end function; 

//same as above, but the output is in the format pattalg

SimpleTriangularAlgebra:=function(n)
	P:=rec<pattalg | index:=[i: i in [1..n]], basis:=[]>;
	for i:=n to 1 by -1 do for k:=i+1 to n do 
		Append(~P`basis,<i,k>);
	end for; end for;
	return P;
end function;


//computes the character degrees of the above algebra
// Uses the "mixed" procedure, starting with pattern algebras. 
// Then treats exceptions.

MTriangularDegrees:=function(n)
	t:=Cputime();
	res,exc,ec:=pattern_fr(SimpleTriangularAlgebra(n),<0,0,0>);
	print "Time=", Cputime(t);
	return AddExceptions(res,exc,ec);
end function;

EmptyAlgebra:=function()
	return rec<algtype | basis:=[], prod:=[], par:=rec<partype| names:=[],eqns:=[],nz:=[],mons:=[]> >;
end function;



// returns the pattern algebra corresponding to the partial order given by C, where C is the list of "adjacent" pairs with respect to the order (i.e. pairs (i,j) such that i<j and no k satisfies i<k<j). 

PatternAlgebra:=function(C)
	repeat
		change:=false;
		for ij in C do for jk in C do
			if ij[2] eq jk[1] and not <ij[1],jk[2]> in C then 
				Append(~C,<ij[1],jk[2]>); change:=true; 
			end if;
		end for; end for;
	until not change;
	A:=EmptyAlgebra();
	A`basis:=C;
	for ij in C do for jk in C do
		if ij[2] eq jk[1] then Append(~A`prod,<ij,jk,<ij[1],jk[2]>,[Integers() | ]>); end if;
	end for; end for;
	return A;
end function;

// same as the previous function, but returns the algebra in the format pattalg

SimplePatternAlgebra:=function(ind,C)
	repeat
		change:=false;
		for ij in C do for jk in C do
			if ij[2] eq jk[1] and not <ij[1],jk[2]> in C then 
				Append(~C,<ij[1],jk[2]>); change:=true; 
			end if;
		end for; end for;
	until not change;
	P:=rec<pattalg | index:=ind, basis:=C>;
	return P;
end function;

AlgebraDegrees:=function(A)
 	return degnonspec(A,<0,0,0>);
end function;

// used by the next function

SortPairs:=function(x,y)
	if x[1] lt y[1] then return 1; 
		elif y[1] lt x[1] then return -1;
		else return x[2]-y[2];
	end if;
end function;

// returns the pattern algebra given by the (totally ordered) set <index> and the relation R obtained by subtracting <zero> from the total order of <index>

PatternFromZeros:=function(index,zero)
	C:=[];
	for i:=1 to #index do for j:=i+1 to #index do if not <index[i],index[j]> in zero then Append(~C,<index[i],index[j]>); end if; end for; end for;
	Sort(~C,SortPairs);
	P:=rec<pattalg | index:=index, basis:=C>;
	return P;
end function;

// Same as MTriangularAlgebra, but works for a general pattern algebra encoded as in the previous function.

DegreesFromZeros:=function(index,zero)
	t:=Cputime();
	res,exc,ec:=pattern_fr(PatternFromZeros(index,zero),<0,0,0>);
	print "Time=", Cputime(t);
	return AddExceptions(res,exc,ec);
end function;

PatternDegrees:=function(ind,C)
	t:=Cputime();
	res,exc,ec:=pattern_fr(SimplePatternAlgebra(ind,C),<0,0,0>);
	print "Time=", Cputime(t);
	return AddExceptions(res,exc,ec);
end function;

