
/* 
 * Magma code associated with "Combinatorial methods of character 
 * enumeration for the unitriangular group."
 *
 * Author: Eric Marberg
 * Version: July 2011.
 * --------------------------------------------------------------
 * The functions below carry out the computations described in 
 * the paper "Combinatorial methods of character enumeration 
 * for the unitriangular group" (http://arxiv.org/abs/1012.2341).  
 * Our functions require the methods defined in Anton Evseev's code,
 * which can be found at
 *
 *   http://math.mit.edu/~emarberg/char_mixed.m
 *
 * This code works in version 2.16-7 of magma, but may not work 
 * in earlier versions.
 *
 * The main methods are
 *
 *   (Constructs the algebras defined in Section 3.1 of the paper)
 *   Get_CrossingAlgebra(Cr);
 *   Get_TildeCrossingAlgebra(Arc,Cr); 
 *
 *   (Attempts to compute polynomials for the nonzero functions N_{Lambda,e}(q))
 *   Compute_N_Lambda_e(Lambda);
 *
 *   (Produces lists of set partitions)
 *   Get_SetPartitions(n, case_n_minus_one);
 *   Get_SetPartitionsUpTo(n);
 *   Get_SetPartitionIterator(n,SP);
 *
 *   (Computes the crossing-connected components of a set partition)
 *   Get_CrConnComponents(Lambda);
 *
 *   (Produces lists of crossing-connected set partitions)
 *   Get_CrConnSetPartitions(n,SP,degree_bound);
 *   Get_CrConnSetPartitionsUpTo(n,SP,degree_bound);
 *
 *   (Produces a map from set partitions to the outputs of Compute_N_Lambda_e)
 *   GetAll_N_Lambda_Map(SP);
 *
 *   (Attempts to compute polynomials for the nonzero functions N_{n,e}(q))
 *   Compute_N_n_e(n,SP,CrConnMap);
 *
 *   (Attempts to compute polynomials for the nonzero functions Tilde_N_{n,e}(q) 
 *    defined in Section 3.4 of the paper. Provides the tables in Appendix A)  
 *   Compute_Tilde_N_n_e(n,e,SP,CrConn);
 *   Compute_Tilde_N_n_e_UpTo(n,e,SP,CrConn);
 *
 *   (Attempts to compute polynomials for the nonzero functions Tilde_N_{n,e}(q) 
 *    defined in Section 3.4 of the paper. Provides the tables in Appendix B)  
 *   Compute_N_n_e_From_Tilde_N_n_e(n_scalar,e,Tilde_N_n_e_Outputs);
 *  
 * See my comments below for detailed information on the usage of these functions.
 * The following example code should also be helpful. If you require further 
 * information, please contact me. 
 *
 */

/*

// ************
// EXAMPLE CODE
// ************

// Make sure to have Evseev's code (http://www.maths.qmul.ac.uk/~evseev/char_mixed.m)
// saved as "char_mixed.m" in the same directory as this file when starting Magma.

load "Combinatorial_Methods_Code.m";


// Define a set partition and compute N_{Lambda,e}(q)
Lambda := {[1,3],[2,4]};
f,fns := Compute_N_Lambda_e(Lambda);

// The output fns is a list of tuples < integer i, polynomial a(x) >
// where N_{Lambda,i}(x) = a(x). In this case, fns is [ <1,x> ];
fns;


// Let's try this with the set partition discussed in Example 3.1 in the paper:
Lambda := { [1,5,7,9,13], [2,6,8,12], [3,10], [4,11] };
f,fns := Compute_N_Lambda_e(Lambda);
fns;

// Here our algorithm fails. However, a slight modification to Evseev's code
// will prevent this failure.  Namely, replace line 742 of Evseev's code: 
//
//   degspec:=function(A,f,N: hist:=[* [], [], [], [] *]) 
// 
// with 
//
//	forward IsAbelian;
//	degspec:=function(A,f,N: hist:=[* [], [], [], [] *])
//
//	if IsAbelian(A) then
//	  succ,res:=proper_count(A`par);
//	  if succ then
//	    return (x-1)^(f[1]+1)*x^f[2]*y^f[3]*res*x^(#A`basis-1), [], []; 
//	  end if; 
//	end if;
//
// The code above should now output in fns the polynomials in Example 3.1.

// NOTE: in all subsequent calculations we assume that Evseev's code does
// not have this modification. We can therefore conclude that the irreducible
// constituents of a supercharacter with shape Lambda are all well-induced
// if calling the method Compute_N_Lambda_e(Lambda) is successful.

// Our algorithm (with Evseev's code unmodified) fails to compute 
// N_{Lambda,e} for only 1 set partition of [12]. Luckily, the
// algorithm does succeed for Lambda^dag where Lambda^dag is the 
// set partition given by applying the map i -> n+1-i to the parts 
// of Lambda, and so we may invoke Proposition 3.2.

// Our algorithm likewise fails to compute N_{Lambda,e} for 34 set 
// partitions of [13].  For 32 of these 34 cases, the algorithm does 
// succeed for Lambda^dag and we may again invoke Proposition 3.2. 
// Of the two remaining cases, one is the set partition 
//
//   { [1,5,7,9,13], [2,6,8,12], [3,10], [4,11] }
//
// described above. The only remaining case is 
//
//   Lambda := {[1,5,9,13],[2,8],[3,7,11],[4,10],[6,12]};
//
// By summing the successful outputs of Compute_N_Lambda_e and subtracting
// from the output of Evseev's function MTriangularDegrees(13); we find
// that N_{Lambda,e}(q) = q^8 - q^7 if e = 17,
//                        q^6 - q^5 + q^4 if e=18,
//                        q^5 - 2*q^4 +2*q^3 - q^2 if e=19,
//                        0 otherwise,
//
// which confirms Theorem 3.3 since each of these functions is a polynomial
// in q-1 with nonnegative integer coefficients.

// One can directly compare the outputs of Compute_N_n_e and Evseev's function 
// MTriangularDegrees with the following code.  

n := 8; // some integer
SP := Get_SetPartitionsUpTo(9);
CrConn := Get_CrConnSetPartitionsUpTo(n,SP,DEGREE_BOUND_OFF);
CrConnMap := GetAll_N_Lambda_Map(CrConn);

f,fns := Compute_N_n_e(n,SP,CrConnMap);
g := MTriangularDegrees(n);

// This difference will be zero for n < 12. For n = 12, 13, one has to
// modify CrConnMap to give the correct outputs for the exceptions
// described above.
f-g;

 
// This set partition is considered in Example 3.3:
Lambda := { [1,4,6], [2,14,15], [3], [5,7,10], [8,11,13], [9], [12] };

// Compute its crossing-connected components
cmpnts := Get_CrConnComponents(Lambda);

f,fns := Compute_N_Lambda_e(Lambda);

g := MV ! 1;
for Gamma in cmpnts do
  st := Standardize(Gamma);
  delta := SetPartition_Degree(Gamma) - SetPartition_Degree(st);
  g := g * Compute_N_Lambda_e(st) * (MV ! y)^delta;
end for;

gns := Get_CountFns(g);

// Now compare fns and gns
fns eq gns;



// The example mentioned in Theorem 3.3(2) is the following:
Lambda := { [1,6,8,13], [2,7,12], [3,9], [4,10], [5,11] };
f,fns := Compute_N_Lambda_e(Lambda);

for pair in fns do
  printf "N_{Lambda,%o}(x+1) = %o\n", pair[1], MV ! Evaluate(pair[2],1,MV ! x + 1);
end for;




// This computes the function Tilde_N_{n,e}(q) defined in Section 3.4.
// It is extremely time consuming to calculate all the functions 
// listed in Appendix A: the bottleneck is the part of the code
// enumerating the set partitions to be considered. As in the next
// example, the best way to deal with this is to implement our
// enumeration algorithms in a faster language and precompute the 
// relevant lists, then load them into Magma.

n := 8;
e := 8;

SP := Get_SetPartitionsUpTo(9);
CrConn := Get_CrConnSetPartitionsUpTo(n+1,SP,e);

// a_vecs is a list of triples <n, i, [a_0,a_1,...,a_ell]> 
// encoding the rows of the tables in Appendix A.
f,fns,a_vecs := Compute_Tilde_N_n_e(n,e,SP,CrConn);

a_vecs;


// The following code enacts the calculations discussed in the
// the proof of Proposition 3.3.  This is an extremely time
// consuming computation for e=6,7,8, due to, among other things, 
// the taxing memory requirements involved in enumerating the 
// crossing-connected set partitions of [n] up to n = 19.
// In practice it is more efficient to precompute these set 
// partitions in a faster language like Java or C (using 
// the same algorithms as those here) and then load the list 
// into Magma.  

e := 4;
SP := Get_SetPartitionsUpTo(9);

for n in [1..2*e+3] do
  CrConn := Get_CrConnSetPartitions(n,SP,e);
  CrConnMap := GetOne_N_Lambda_Map(CrConn);

  printf "*** RESULT: All crossing-connected set partitions Lambda of [%o] have N_{Lambda,e}(q) = 0 when e < %o ***\n\n",
         n,N_Lambda_Map_MinimumDegree(CrConnMap);
end for;


// This computes the rows in the tables in Appendix B.
// As in the previous code, this is an extremely time
// consuming calculation for e = 6,7,8.

e := 4; 

SP := Get_SetPartitionsUpTo(9);
CrConn := Get_CrConnSetPartitionsUpTo(2*e+2,SP,e);
Tilde_N_n_e_Outputs := Compute_Tilde_N_n_e_UpTo(2*e+1,e,SP,CrConn);

// a_vecs is a list of triples <e, i, [a_0,a_1,...,a_ell]> 
// encoding the rows of the tables in Appendix B.
// f is N_{10,e}(q).
f,fns,a_vecs := Compute_N_n_e_From_Tilde_N_n_e(10,e,Tilde_N_n_e_Outputs);

a_vecs;
f;

// *******************
// END OF EXAMPLE CODE
// *******************

*/







load "char_mixed.m";

/*
 * Function: IsAbelian
 * ---------------------------------------------------------
 * Returns true if the input algebra A is abelian, and false 
 * otherwise.
 *
 */

IsAbelian:=function(A) 

  for ijkl in A`prod do
    if not <ijkl[2], ijkl[1], ijkl[3], ijkl[4]> in A`prod then
      return false;
    end if;
  end for;

  return true;

end function;


/*
 * Function: Get_ArcSet
 * ----------------------------------------------------------------
 * Given a set partition Lambda, this function returns Arc(Lambda).
 * The set partition Lambda must be represented as a list of lists 
 * of numbers in increasing order. For example,
 *
 *   Lambda := {[1,4,5],[2,3],[6],[7,9],[8]}; 
 *   Get_ArcSet(Lambda);
 *
 * returns the list [<1,4>, <4,5>, <2,3>, <7,9>].
 *
 */

Get_ArcSet:=function(Lambda)

  Arc := [];
  for B in Lambda do

    for i in [1..#B-1] do
      Append(~Arc,<B[i],B[i+1]>);
    end for;

  end for;

  return Arc;

end function;


/*
 * Function: Get_CrossingSet
 * --------------------------------------------------------------
 * Given Arc(Lambda) for a set partition Lambda, this function 
 * returns Cr(Lambda), the set of pairs (i,j) for which there are 
 * i<j<k<l with (i,k),(j,l) in Arc(Lambda). For example,
 *
 *   Lambda := {[1,4,5],[2,3],[6],[7,9],[8]}; 
 *   Arc := Get_ArcSet(Lambda);
 *   Get_CrossingSet(Arc);
 *
 * returns the list [<1,4>, <4,5>, <2,3>, <7,9>].
 *
 */

Get_CrossingSet:=function(Arc)

  Cr := [];
  for a in Arc do
    for b in Arc do

      i := a[1];
      j := b[1];
      k := a[2];
      l := b[2];

      if i lt j and j lt k and k lt l then
        Append(~Cr,<i,j>);
      end if;

    end for; 
  end for;

  return Cr;

end function;


/*
 * Function: Get_CrossingAlgebra
 * --------------------------------------------------------------
 * This function returns an algebra over the integers which is
 * isomorphic to the quotient of a pattern algebra by an ideal
 * which is a pattern algebra.
 *
 * (Recall, a pattern algebra is an algebra spanned by elementary
 * matrices.)
 *
 * The input C is a list of pairs <i, j> representing positions 
 * in a matrix above the diagonal.  The set of these pairs must 
 * be of the form J - K for two sets of J supset K, which have 
 * the property that the matrices e_ij for (i,j) in J span a 
 * pattern algebra A and the matrices e_ij for (i,j) in K span an 
 * ideal of A.  If C is not of this form then the output algebra
 * may fail to be associative.
 *
 */

Get_CrossingAlgebra:=function(C)

  A:=EmptyAlgebra();
	

  // max = the maximum value of j over all <i,j> in C.

  max := 0;
  for ij in C do
    if ij[2] gt max then max := ij[2]; 
    end if;
  end for;


  // The basis {b_i} of A just consists of the elements of C,
  // but to work with Evseev's functions these must be listed 
  // in an order such that each product b_i*b_j is a linear 
  // combination of b_k's with k > max(i,j). To accomplish this,
  // we add the elements of C to A`basis by iterating first by 
  // diagonal d then by row i.

  for d in [1..max-1] do
    for i in [1..max-d] do

      if <i,i+d> in C then
        Append(~A`basis,<i,i+d>);
      end if;

    end for; 
  end for;


  // Now we specify the product in A

  for kl in A`basis do 
    for lm in A`basis do

      if kl[2] eq lm[1] and <kl[1],lm[2]> in C then
        Append(~A`prod,<kl,lm,<kl[1],lm[2]>,[Integers()| ]>);
      end if;

    end for; 
  end for;

  return A;

end function;


/*
 * Function: Get_TildeCrossingAlgebra
 * --------------------------------------------------------
 * Given the arc set Arc = Arc(Lambda) and the crossing set 
 * Cr = Cr(Lambda) of a set partition Lambda of [n], this 
 * function returns the algebra over Z with the basis 
 * 
 *   { e_ij : (i,j) in Cr } cup { z_Lambda }
 *
 * and the multiplication defined by equation (3.1) in 
 * "Combinatorial methods ...". Tensoring this algebra with 
 * the finite field FF_q gives the algebra
 *
 *   $\widetilde \mathfrak{C}_\Lambda(q)$ 
 *
 * defined in Section 3.1 of the paper.
 *
 */

Get_TildeCrossingAlgebra:=function(Arc,Cr)

  // The output algebra is formed by adding 
  // one additional central basis element to
  // to Get_CrossingAlgebra(Cr).

  z := Arc[1];
  A:=Get_CrossingAlgebra(Cr);
  Append(~A`basis,z);
	
  for kl in A`basis do 
    for lm in A`basis do
      if kl[2] eq lm[1] and <kl[1],lm[2]> in Arc then
        Append(~A`prod,<kl,lm,z,[Integers()| ]>);
      end if;
    end for; 
  end for;

  return A;

end function;


/*
 * Function: SetPartition_Degree, SetPartition_MinConstitDegree
 * --------------------------------------------------------------
 * Given a set partition Lambda, the function SetPartition_Degree 
 * returns the integer d(Lambda) defined in Section 2.5 of 
 * "Combinatorial methods...". Recall that any supercharacter of 
 * UT_n(q) with shape Lambda has degree q^d(Lambda).
 *
 * Given a set partition Lambda, SetPartition_MinConstitDegree
 * returns the integer d(Lambda) - |Cr(Lambda)|. If e is less
 * than this integer then N_{Lambda,e}(q) = 0.
 *
 */

SetPartition_Degree:=function(Lambda)

  d := 0;

  for part in Lambda do
    for i in [1..(#part - 1)] do
      d := d + part[i+1] - part[i] - 1;
    end for;
  end for;

  return d;

end function;


SetPartition_MinConstitDegree:=function(Lambda)
  return SetPartition_Degree(Lambda) - #Get_CrossingSet(Get_ArcSet(Lambda));
end function;


/*
 * Function: Get_CountFns
 * -------------------------------------------------------------------
 * Given a bivariate polynomial F in x and y, this method returns the 
 * list of all pairs <e, f(x)> where e is a nonnegative integer such  
 * that the coefficient of y^e in F(x,y) is nonzero and equal to f(x).
 *
 */

Get_CountFns:=function(f);

  count_fns := [];
  for i in [0..MV!Degree(f,2)] do
    g := MV ! Coefficient(f,2,i);
    if g ne 0 then Append(~count_fns,<i,g>);
    end if;
  end for;

  return count_fns;

end function;


/*
 * Function: Compute_N_Lambda_e
 * ------------------------------------------------------------------
 * Given a set partition of [n], represented as a list of lists of
 * integers in increasing order (e.g., [ [1,5,6], [2,4], [3] ]), this
 * function attempts to compute polynomials which are equal to the 
 * functions N_{Lambda,e}(x) defined in "Combinatorial methods ..." 
 * for all integers e.  
 *
 * If the algorithm is successful:
 *  -The first output is a polynomial in x and y, in which the 
 *   coefficient of y^e is N_{Lambda,e}(x).
 *  -The second output is the list of all pairs <e, f(x)> where e is 
 *   a nonnegative integer such that N_{Lambda,e}(x) = f(x) != 0.  
 *
 * If the algorithm fails, the outputs are 0, [].
 *
 */

Compute_N_Lambda_e:=function(Lambda)

  Arc:=Get_ArcSet(Lambda);
  Cr:=Get_CrossingSet(Arc);
  d_Lambda := SetPartition_Degree(Lambda);

  
  // If #Cr is zero then any supercharacter with shape Lambda
  // is irreducible, so N_{Lambda,e}(q) = 1 if e = d(Lambda)
  // and is zero for all other values of e.

  if #Cr eq 0 then
    f := (MV ! y)^d_Lambda;
    return f,Get_CountFns(f);
  end if;
  

  // A supset B are algebras over the integers Z.
  //
  // With respect to the notation in "Combinatorial methods ...":
  //   A otimes FF_q is isomorphic to $\widetilde\mathfrak{C}_\Lambda(q)$,
  //   B otimes FF-q is isomorphic to $\mathfrak{C}_\Lambda(q)$.

  A := Get_TildeCrossingAlgebra(Arc,Cr);
  B := Get_CrossingAlgebra(Cr);


  // Evseev's function AlgebraDegrees takes as input a nilpotent 
  // algebra N over Z of a certain special form; our algebras A and
  // B are defined to be of this form.  AlgebraDegrees attempts
  // to compute a bivariate polynomial in x and y, such that the
  // the coefficient of y^e is a polynomial f(x) whose value at 
  // prime powers x = q is the number of irreducible characters 
  // of degree q^e of the algebra group 1 + N otimes FF_q.
  //
  // AlgebraDegrees has three outputs: the first is this polynomial,
  // and the last two record data associated with the characters 
  // which the algorithm failed to count. In particular, Evseev's
  // algorithm is successful if and only if both of the second two
  // outputs are empty.

  res1,failed_a_1,failed_b_1 := AlgebraDegrees(A);
  res2,failed_a_2,failed_b_2 := AlgebraDegrees(B);

  if #failed_a_1 ne 0 or #failed_b_1 ne 0 or #failed_a_2 ne 0 or #failed_b_2 ne 0 then 

    printf "Failed to compute N_{Lambda,e}(q) for Lambda =\n";
    print Lambda; print "";

    return 0,[];
  end if;
 

  // Assuming Evseev's algorithm was successful in both cases, 
  // Theorem 3.1 in "Combinatorial methods ..." asserts that 
  //
  //   N_{Lambda,e}(x) = ( a(x) - b(x) ) / ( x - 1 ) 
  //
  // where a(x) is the coefficient of y^f in res1 and b(x) is 
  // the coefficient of y^f in res2, for f = #Cr - d_Lambda + e.
  // This accounts for our return value.

  f := (res1-res2) div (MV ! x - 1) * (MV ! y)^(d_Lambda-#Cr);
  return f,Get_CountFns(f);

end function;


/*
 * Function: Get_SetPartitions, Get_SetPartitionExtensions
 * -----------------------------------------------------------------
 * A set partition is a set of disjoint lists of integers arranged
 * in increasing order. E.g., { [1,3,5,6], [2,4,8], [7] }.
 *
 * Given n and the list of all set partitions of {1,2,...,n-1}, the
 * function Get_SetPartitions  outputs the list of all set partitions 
 * of [n] = {1,2,..,n}. For example,
 *
 *   SP_0 := [{}];
 *   SP_1 := Get_SetPartitions(1,SP_0);
 *   SP_2 := Get_SetPartitions(2,SP_1);
 *
 * and so on.
 *
 * The function Get_SetPartitionExtensions is a helper method which
 * takes as input an integer n and set partition Lambda of [n-1],
 * and which returns the list of set partitions of [n] formed by 
 * inserting n somewhere in Lambda.
 *
 */

Get_SetPartitionExtensions:=function(n,Lambda);
  extns := [Include(Lambda,[n])];
  for part in Lambda do
    Mu := Lambda;
    Exclude(~Mu,part);
    Include(~Mu,Append(part,n));
    Append(~extns, Mu);
  end for;
  return extns;
end function;

Get_SetPartitions:=function(n, case_n_minus_one);

  // Each set partition of [n] is obtained from a unique set partition
  // of [n-1], either by placing n in its own singleton part, or by
  // adding n to one of the existing parts.

  case_n := [];
  for Lambda in case_n_minus_one do
    case_n := case_n cat Get_SetPartitionExtensions(n,Lambda);
  end for;

  return case_n;

end function;


/*
 * Function: Get_SetPartitionsUpTo
 * --------------------------------------------------------
 * Given a positive integer n, this method returns the list
 * [SP_1, SP_2, ..., SP_n] where 
 *
 *   SP_0 := [{}];
 *   SP_1 := Get_SetPartitions(1,SP_0);
 *   SP_2 := Get_SetPartitions(2,SP_1);
 *
 * and so on.
 *
 */

Get_SetPartitionsUpTo:=function(n);

  printf "\n";
  SP := [];

  prev := [{}];
  for i in [1..n] do
    printf "Computing set partitions of [%o] . . . ", i;
    prev := Get_SetPartitions(i,prev);
    Append(~SP,prev);
    printf "# = %o\n", #prev;
  end for;

  printf "\n";
  return SP;

end function;


/*
 * Function: Standardize
 * ----------------------------------------------------
 * Given a set partition Lambda, this method return its
 * standardization.  For example,
 *
 *   Lambda := {[2,5,6], [3,9], [11]};
 *   Standardize(Lambda);
 *
 * returns {[1,3,4], [2,5], [6]};
 *
 */

Standardize:=function(Lambda);

  numbers := [];
  for part in Lambda do
    numbers := numbers cat part;
  end for;

  Sort(~numbers);
  S := [];
  for i in [1..#numbers] do
    Append(~S,<numbers[i],i>);
  end for;

  m := map< numbers -> [1..#numbers] | S>;

  Gamma := {};
  for part in Lambda do

    new_part := [];
    for i in part do 
      Append(~new_part,i@m);
    end for;

    Include(~Gamma,new_part);

  end for;

  return Gamma;

end function;


/*
 * Function: Get_CrConnComponents, InsertArc
 * ------------------------------------------------------------
 * Given a set partition Lambda as input, this method returns a 
 * list of the partitions crossing-connected components, as 
 * defined in "Combinatorial methods ...". 
 *
 * InsertArc is a helper function used in Get_CrConnComponents.
 *
 */

InsertArc:=function(a,comp);

  left_join := [];
  right_join := [];
  updated_comp := {};

  for part in comp do

    // Case: comp has a part of of the form [a[2], x, y, z, ...]
    if part[1] eq a[2] then 
      right_join := Remove(part,1);

    // Case: comp has a part of of the form [..., x, y, z, a[1]]
    elif part[#part] eq a[1] then
      left_join := Remove(part,#part);

    else 
      Include(~updated_comp,part);
    end if;

  end for;

  Include(~updated_comp,left_join cat [a[1],a[2]] cat right_join);
  return updated_comp;

end function;


Get_CrConnComponents:=function(Lambda);

  components := [];

  // Add all singleton components
  for part in Lambda do
    if #part eq 1 then Append(~components,{part});
    end if;
  end for;

  // For each arc in Arc(Lambda), we determine its equivalence class
  // of arcs generated by the relation (i,j) ~ (k,l) if (i,j) crosses
  // (k,l), and then form the crossing-connected component of Lambda 
  // corresponding to this set of arcs.  
  Arc := Get_ArcSet(Lambda); 
  while #Arc ne 0 do

    next_comp := {};
    to_consider := [Arc[1]];    
    Remove(~Arc,1);
   
    while #to_consider ne 0 do
    
      // This is the set of arcs crossing the current arc.
      // On the next iteration we must compute the arcs 
      // which cross each of these.
      to_consider_next := [];

      for a in to_consider do 

        // Insert the arc a into component under construction
        next_comp := InsertArc(a, next_comp);

        // We iterate from the end of Arc to its beginning so that
        // we can remove elements of Arc as we go along.
        len := #Arc;
        for i in [0..(len-1)] do

          index := len - i;
          b := Arc[index];

          // Check if a crosses b, and if so add this to the list 
          // of arcs to be considered on the next iteration
          if (b[1] lt a[1] and a[1] lt b[2] and b[2] lt a[2]) or
             (a[1] lt b[1] and b[1] lt a[2] and a[2] lt b[2]) then
            Append(~to_consider_next,b) ;
            Remove(~Arc,index);
          end if;

        end for;
      end for;

      to_consider := to_consider_next;

    end while;

    Append(~components,next_comp);

  end while;

  return components;

end function;


/*
 * Function: Get_SetPartitionIterator, Get_CrConnSetPartitionIterator
 * ------------------------------------------------------------------
 * The following functions are used to deal with memory limitations
 * which arise from having to work with very large lists of set 
 * partitions.
 *
 * Get_SetPartitionIterator takes as input 
 *   n = a positive integer, 
 *   SP = the return value of Get_SetPartitionsUpTo(m); where m is 
 *        any positive integer,
 * 
 * The method returns three outputs: it, PARAMS, VALUE. The first output
 * is a procedure, which is used as follows:
 *
 *   it, PARAMS, VALUE := Get_SetPartitionsIterator(n,SP);
 *   while true do
 *   
 *     it(~PARAMS,~VALUE);
 *     if VALUE eq SetPartitionIterator_DONE_FLAG then break;
 *     end if;
 *
 *     ...
 *
 *   end while;
 *
 * The while loop in this example iterates over all set partitions of [n],
 * but only requires the set partitions of [m] for some m <= n be stored 
 * in memory.
 *
 * Get_CrConnSetPartitionIterator takes as input 
 *   n = a positive integer, 
 *   SP = the return value of Get_SetPartitionsUpTo(m); where m is 
 *        any positive integer,
 *   degree_bound = a nonnegative integer
 *
 * Its outputs and usage are the same as for Get_SetPartitionIterator, 
 * except that only set partitions Lambda with d(Lambda) - |Cr(Lambda)| 
 * <= degree_bound are assigned to VALUE. When DEGREE_BOUND_OFF is the 
 * inputted, the function treats degree_bound as positive infinity. 
 * This function is used exclusively in the method Get_CrConnSetPartitions.
 *
 * Get_SetPartitionIterator_Helper is a helper function for these methods.
 *
 */

DEGREE_BOUND_OFF := -1;
SetPartitionIterator_DONE_FLAG := {};

Get_SetPartitionIterator_Helper:=function(n,SP,restrictive_flag,cr_conn_flag);

  m := Min(#SP,n);
  SP := SP[m];

  it:=procedure(~PARAMS,~VALUE);

    while PARAMS[1] eq #PARAMS[3] do

      PARAMS[1] := 0;
      PARAMS[2] := PARAMS[2] + 1;

      if PARAMS[2] gt # SP then
        VALUE := SetPartitionIterator_DONE_FLAG;
        break;
      end if;

      PARAMS[3] := [SP[PARAMS[2]]];

      for i in [m+1..n] do
      
        new_block := [];
        for Lambda in PARAMS[3] do

          to_add := true;

          if restrictive_flag ne DEGREE_BOUND_OFF then

            bound := SetPartition_MinConstitDegree(Lambda);

            // We increase the bound when cr_conn_flag is true because
            // in that case the output set partitions will be shifted 
            // to the right by one column, which increases the value
            // of d(Lambda) - |Cr(Lambda)|

            if cr_conn_flag then 
              bound := bound + #Get_ArcSet(Lambda) - (i-1 -#Lambda);
            end if;

            to_add := bound le restrictive_flag;

          end if;

          if to_add then
            new_block := new_block cat Get_SetPartitionExtensions(i,Lambda);
          end if;

        end for;

        PARAMS[3] := new_block;

      end for;

    end while;

    PARAMS[1] := PARAMS[1] + 1;

    if VALUE ne SetPartitionIterator_DONE_FLAG then 
      VALUE := PARAMS[3][PARAMS[1]];
    end if;

  end procedure;

  return it, [* 0,0,[] *], {1};

end function;


Get_SetPartitionIterator:=function(n,SP);
  return Get_SetPartitionIterator_Helper(n,SP,DEGREE_BOUND_OFF,false);
end function;


Get_CrConnSetPartitionIterator:=function(n,SP,degree_bound);
  return Get_SetPartitionIterator_Helper(n,SP,degree_bound,true);
end function;


/*
 * Function: ArcSet_to_SetPartition
 * ------------------------------------------------------------
 * Given a list arcs of pairs <i,j> and integers lb <= ub, this
 * function returns the set partition of [lb..ub] whose arc set 
 * is arcs.
 *
 */

ArcSet_to_SetPartition:=function(arcs,lb,ub)

  Gamma := {};
  for i in [lb..ub] do
    Include(~Gamma,[i]);
  end for;

  for a in arcs do
    Gamma := InsertArc(a,Gamma);
  end for;

  return Gamma;

end function;


/*
 * Function: SetPartition_Shift
 * -----------------------------------------------------------
 * This function "shifts" the matrix of an input set partition
 * Lambda. The others inputs lb, ub, dx, dy are integers. What
 * is returned is the unique set partition of {lb,lb+1,...,ub}
 * whose arc set is { (i+dx, j+dy) : (i,j) in Arc(Lambda) }.
 *
 */

SetPartition_Shift:=function(lb,ub,dx,dy,Lambda)

  Arc := Get_ArcSet(Lambda);
  new_arc := [];

  for a in Arc do
    Append(~new_arc, <a[1]+dx,a[2]+dy>);
  end for;

  return ArcSet_to_SetPartition(new_arc,lb,ub);

end function;


/*
 * Function: PrintProcedure
 * ---------------------------------------------------------------
 * Helper procedure used to print tracking output in a few places.
 *
 */

PrintProcedure:=procedure(~dot,incr,total);

  if dot mod incr eq 0 then printf ". "; end if;
  dot := dot + 1;
  if dot mod (32*incr) eq 0 then print Floor(100*dot / total), "%"; end if;

end procedure;


/*
 * Function: Get_CrConnSetPartitions, Get_CrConnSetPartitions_Helper, ProcessLambda
 * --------------------------------------------------------------------------------
 * The function Get_CrConnSetPartitions takes as input a positive integer n, 
 * a list whose ith entry is the list of all set partitions of [i], and a 
 * nonnegative integer degree_bound.  
 * 
 * It returns a list of all crossing-connected set partitions Lambda of [n] 
 * with d(Lambda) - |Cr(Lambda)| <= degree_bound.  (If this inequality fails 
 * then N_{Lambda,e}(q) = 0 for all e <= degree_bound.)  If DEGREE_BOUNDED_OFF
 * is inputted for degree_bound, then the method treats degree_bound as positive 
 * infinity.
 *
 * ProcessLambda is a helper function used to construct all crossing-connected 
 * set partitions which can be formed by adding 1 and n to a set partition of 
 * [2..n-1].
 *
 * Get_CrConnSetPartitions_Helper is the same function as Get_CrConnSetPartitions
 * except that it takes an additional input to determine whether it prints any
 * status updates.
 *
 */

ProcessLambda:=procedure(Lambda,n,degree_bound,~cr_conn_set);

  Lambda := SetPartition_Shift(2,n-1,1,2,Lambda);

  singletons := 0;
  for part in Lambda do
    if #part eq 1 then 
      singletons := singletons + 1;
    end if;
  end for;

  // We are adding 1 and n to Lambda is all possible ways. If 
  // Lambda has more than two parts with a single element, then
  // we will be left with at least one singleton part, in which
  // case the result is not crossing-connected.

  if singletons le 2 then
    for part_1 in Lambda do

      Gamma := Lambda;
      Exclude(~Gamma, part_1);
      Include(~Gamma,[1] cat part_1); 
 
      if (degree_bound eq DEGREE_BOUND_OFF) or (SetPartition_MinConstitDegree(Gamma) le degree_bound) then

        for part_n in Gamma do

          Kappa := Gamma;
          Exclude(~Kappa, part_n);
          Include(~Kappa,part_n cat [n]);
     
          if  #Get_CrConnComponents(Kappa) eq 1 then
            Append(~cr_conn_set,Kappa);
          end if;

        end for;
      end if;
    end for; 
  end if;

end procedure;


Get_CrConnSetPartitions_Helper:=function(n,SP,degree_bound,OUTPUT_ON);

  t:=Cputime();	
  cr_conn_set := [];

  // If n is small enough then we don't both with complicated extensions
  // and just produce the crossing-connected partitions from a brute 
  // force search.

  if n lt 4 then 
    for Lambda in SP[n] do
      if #Get_CrConnComponents(Lambda) eq 1 then
        Append(~cr_conn_set,Lambda);
      end if;
    end for;
    print "";
    return cr_conn_set;
  end if; 

  incr := 20;
  total := #SP[Min(#SP,n-3)];

  if OUTPUT_ON then printf "\nEach '.\' represents %o iterations of %o\n", incr, total;
  end if;

  // We input n-3 to the iterator because crossing-connected set partitions 
  // never have i,i+1 in the same part, and because they must have an arc
  // of the form <1,i> and <j,n> for some i,j.  Therefore we can shift 
  // a set partition of n-3 to become a set partition of [2..n-1] and then
  // examine all the ways to add arcs involving 1 and n to find the 
  // crossing-connected partitions of interest.
 
  it,PARAMS,Lambda := Get_CrConnSetPartitionIterator(n-3,SP,degree_bound);
  dot := 0;

  while true do

    prev_t_index := PARAMS[2];
    it(~PARAMS,~Lambda);

    if OUTPUT_ON and prev_t_index lt PARAMS[2] then PrintProcedure(~dot,incr,total);
    end if;

    if Lambda eq SetPartitionIterator_DONE_FLAG then break; end if;
    ProcessLambda(Lambda,n,degree_bound,~cr_conn_set);
    
  end while;

  if OUTPUT_ON then print "\nTime=", Cputime(t),"\n";
  end if;

  return cr_conn_set;

end function;


Get_CrConnSetPartitions:=function(n,SP,degree_bound);
  return Get_CrConnSetPartitions_Helper(n,SP,degree_bound,true);
end function;


/*
 * Function: Get_CrConnSetPartitionsUpTo
 * ----------------------------------------------------------------
 * Returns the list [CrConn_1,CrConn_2,...,CrConn_n] where CrConn_i 
 * is the output of Get_CrConnSetPartitions(i,SP,degree_bound).
 *
 */

Get_CrConnSetPartitionsUpTo:=function(n,SP,degree_bound);

  printf "\n";
  CrConn := [];

  if degree_bound gt DEGREE_BOUND_OFF then
    printf "*** Irreducible constituent degree bound = q^%o ***\n\n",degree_bound;
  else
    printf "*** Irreducible constituent degree bound = infinity ***\n\n";
  end if;

  for i in [1..n] do
    printf "Computing crossing-connected set partitions of [%o] . . . ", i;
    Append(~CrConn,Get_CrConnSetPartitions_Helper(i,SP,degree_bound,true));
  end for;

  printf "\n";
  return CrConn;

end function;


/*
 * Function: GetOne_N_Lambda_Map, GetAll_N_Lambda_Map
 * ---------------------------------------------------------------
 * Given a list of set partitions, GetOne_N_Lambda_Map returns the
 * the map assigning each set partition Lambda to the first output
 * of Get_N_Lambda_e(Lambda). 
 *
 * If SP = [SP1,SP2,SP3,...] is a list of pairwise disjoint lists 
 * of set partitions, then GetAll_N_Lambda_Map(SP) returns the 
 * output of GetOne_N_Lambda_Map(SP1 cat SP2 cat SP3 cat ...).
 *
 */

GetOne_N_Lambda_Map:=function(set);

  dom := [];
  codom := [];
  S := [];

  for Lambda in set do
    f := Compute_N_Lambda_e(Lambda);
    Append(~dom,Lambda);
    Append(~codom,f);
    Append(~S,<Lambda,f>);
  end for;

  return map< dom -> codom | S>, S;

end function;


GetAll_N_Lambda_Map:=function(SP);

  dom := [];
  codom := [];
  S := [];

  for i in [1..#SP] do

    set := SP[i];
    printf "Computing N_Lambda_Map for input list, iteration %o . . . ", i;

    m,S_tmp := GetOne_N_Lambda_Map(set);
    dom := dom cat Domain(m);
    codom := codom cat Codomain(m);
    S := S cat S_tmp;

    printf "#domain = %o\n", #S;

  end for;

  return map< dom -> codom | S>;

end function;


/*
 * Function: N_Lambda_Map_MinimumDegree
 * -----------------------------------------------------------------
 * Let m be a map assigning some set of set partitions Lambda to the
 * polynomials given as the first output of Get_N_Lambda_e(Lambda).
 *
 * This function iterates over all set partitions Lambda in the 
 * domain of m and returns the least integer e such that some
 * Lambda has N_{Lambda,e}(q) neq 0.
 *
 * If the domain of m is the set of crossing-connected set 
 * set partitions of [n], then the output of this function is:
 *
 *       n = 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19
 *  Output = 0, 0, 0, 1, 2, 2, 3, 3, 4, 4,  5,  5,  6,  6,  7,  7,  8,  8,  9
 *
 * This is an extremely time consuming calculation. In practice,
 * we implemented the function Get_CrConnSetPartitions in Java to
 * first enumerate the crossing-connected set partitions more
 * rapidly, and then used our functions here to compute the 
 * corresponding maps.
 *
 */

N_Lambda_Map_MinimumDegree:=function(m);

  ans := MV ! 0;
  for elem in Codomain(m) do
    ans := ans + elem;
  end for;

  for i in [1..MV ! Degree(ans,2)] do
    if MV ! Coefficient(ans,2,i) ne 0 then
      return i;
    end if;
  end for;

  return 0;

end function;


/*
 * Function: Compute_N_Lambda_e_From_Map
 * -----------------------------------------------------------------
 * The inputs to this function are Lambda, a set partition, and map,
 * a map whose domain is a set of crossing-connected set partitions,
 * and whose values are the first outputs of Compute_N_Lambda_e.
 *
 * This function computes the first output of Compute_N_Lambda_e
 * for the into Lambda, by decomposing Lambda into its crossing-connected
 * components, looking up their values in map, and then invoking 
 * Theorem 3.2 in "Combinatorial methods ...".
 *
 * If any of the crossing-connected components of Lambda do not appear
 * in the domain of map, the function returns 0.
 */

Compute_N_Lambda_e_From_Map:=function(Lambda,map)

  next := MV ! 1;
  cmpnts := Get_CrConnComponents(Lambda);
    
  for Gamma in cmpnts do
    st := Standardize(Gamma);

    if not (st in Domain(map)) then return 0;
    end if;

    delta := SetPartition_Degree(Gamma) - SetPartition_Degree(st);
    next := next * st@map * (MV ! y)^delta;
  end for;

  return next;

end function;


/*
 * Function: Compute_N_n_e
 * ----------------------------------------------------------------
 * This function attempts to compute polynomials which are equal to 
 * the functions N_{n,e}(x) giving the number of irreducible 
 * characters of degree x^e of the n-by-n unitriangular group over 
 * a finite field with x elements.
 *
 * The inputs are:
 *   n = a positive integer
 *
 *   SP = the output from Get_SetPartitionsUpTo(m); where m is ANY
 *        integers between 0 and n.
 *
 *   CrConnMap = the output from GetAll_N_Lambda_Map(CrConn); where
 *               CrConn is the output from the function
 *               Get_CrConnSetPartitionsUpTo(n,SP,DEGREE_BOUND_OFF);
 *
 * If successful, the outputs are:
 *  -A polynomial in x and y, with N_{n,e}(x) = coefficient of y^e.
 *  -The list of all pairs <e, f(x)> where e is a nonnegative 
 *   integer such that N_{n,e}(x) = f(x) != 0.  
 *
 * If the algorithm fails, the outputs are 0, [].
 *
 * NOTE: This computes the same output as Evseev's function 
 *       MTriangularDegrees(n); and is significantly less efficient
 *       (about ten times slower).
 */

bell := [1, 2, 5, 15, 52, 203, 877, 4140, 21147, 115975, 678570, 4213597, 27644437, 
         190899322, 1382958545, 10480142147, 82864869804, 682076806159, 5832742205057,  
         51724158235372];

Compute_N_n_e:=function(n,SP,CrConnMap);

  t:=Cputime();	    

  incr := 1000;
  total := bell[n];

  printf "\nEach '.\' represents %o iterations of %o\n", incr, total;

  dot := 0;

  // Iterates over all set partitions Lambda of [n] and 
  // sums the resulting functions N_{Lambda,e}(q)

  it,PARAMS,Lambda := Get_SetPartitionIterator(n,SP);

  f := MV ! 0;
  while true do

    it(~PARAMS,~Lambda);
    if Lambda eq SetPartitionIterator_DONE_FLAG then break;
    end if;

    // Computes crossing-connected components of Lambda, then invokes 
    // Theorem 3.2 in "Combinatorial methods ..." to compute the functions
    // N_{Lambda,e}(q) from the stored values in CrConnMap.

    next := Compute_N_Lambda_e_From_Map(Lambda,CrConnMap);

    if next eq 0 then return 0,[]; end if;
    f := f + (MV ! x - 1)^(n - #Lambda) * next;

    // Some output print commands
    PrintProcedure(~dot,incr,total);

  end while;

  print "\nTime=", Cputime(t),"\n";
  return f, Get_CountFns(f);

end function;


/*
 * Function: Get_Tilde_N_n_e_Iterator, Tilde_N_n_e_Iterator_Setup, SetupCompositions
 * ---------------------------------------------------------------------------------
 * Get_Tilde_N_n_e_Iterator takes four inputs:
 *   Lambda = a crossing-connected set partition of [n]
 *   weak_comp = a weak composition with n-1 parts
 *   SP = the output from Get_SetPartitionsUpTo(m); for some m
 *   degree_bound = a positive integer; the set partitions A iterated over 
 *                have d(A) -|Cr(A)| <= degree_bound
 *
 * Let S be the set of integers { i + sum_{j<i} weak_comp[j] : i in [1..n] }.
 * From the input crossing-connected set partition Lambda and the weak composition
 * one obtains a crossing-connected set partition Gamma of S, by applying the order
 * preserving bijection [1..n] -> S to Lambda.
 *
 * The function returns an iterator in the sense of Get_SetPartitionIterator above
 * which iterates over all set partitions of [ n+|weak_comp| ] which have 
 * Gamma as a crossing-connected component.
 *
 * Tilde_N_n_e_Iterator_Setup and SetupCompositions are helper methods.
 *
 */

SetupCompositions:=function(N,full_frame)

  comp := [];
  comp_shifts := [];

  marks := [];
  for i in [1..N] do
    Append(~marks,"X");
  end for;

  for part in full_frame do
    if #part eq 1 then 
      marks[part[1]] := "S";
    else 
      marks[part[1]] := "L";
      marks[part[#part]] := "R";
    end if;
  end for;


  i := 1;
  while i le N do
  
    if marks[i] eq "S" or marks[i] eq "R" then
      
      next_comp_part := 1;
      Append(~comp_shifts,i-1);
      
      while i lt N do
        i := i + 1;

        if marks[i] eq "R" or marks[i] eq "X" then 
          i := i - 1;
          break; 
        elif marks[i] eq "L" then 
          next_comp_part := next_comp_part + 1;
          break;
        else next_comp_part := next_comp_part + 1;
        end if;
      end while;

      Append(~comp,next_comp_part);
 
    end if;

    i := i + 1;

  end while;

  return comp,comp_shifts;

end function;


Tilde_N_n_e_Iterator_Setup:=function(Lambda,weak_comp,SP,degree_bound)

  n := 0;
  for part in Lambda do n := n + #part; end for;

  frame_shifts := [1];
  sum := 0;

  for i in [2..n] do
    sum := sum + weak_comp[i-1];
    Append(~frame_shifts,i+sum);
  end for;

  frame := {};
  for part in Lambda do
    new_part := [];
    for i in part do
      Append(~new_part, frame_shifts[i]);
    end for;
    Include(~frame,new_part);
  end for;

  full_frame := ArcSet_to_SetPartition(Get_ArcSet(frame),1,n+sum);
  comp,comp_shifts := SetupCompositions(n+sum,full_frame);


  iter := [* *];
  param := [* *]; 
  value := [* *];

  for i in [1..#comp] do
    nxt_i, nxt_p, nxt_v := Get_SetPartitionIterator_Helper(comp[i],SP,degree_bound,false);
   
    if i eq 1 then Append(~value,nxt_v);
    else nxt_i(~nxt_p,~nxt_v); Append(~value,nxt_v); 
    end if;

    Append(~iter,nxt_i);
    Append(~param,nxt_p);

  end for;


  return n+sum, comp, comp_shifts, frame, full_frame, [*iter,param,value*];

end function;


Tilde_N_n_e_Iterator_DONE_FLAG := -1;

Get_Tilde_N_n_e_Iterator:=function(Lambda,weak_comp,SP,degree_bound)

  N, comp, comp_shifts, frame, full_frame, rvals := Tilde_N_n_e_Iterator_Setup(Lambda,weak_comp,SP,degree_bound);
  
  n_Lambda_e,fns := Compute_N_Lambda_e(full_frame);
  skip_bool := #fns gt 0 and fns[1][1] gt degree_bound;

  it := procedure(~PARAMS,~VALUE)

    if skip_bool then VALUE := Tilde_N_n_e_Iterator_DONE_FLAG;
    else 

      iter := PARAMS[1];
      param := PARAMS[2];
      value := PARAMS[3];

      for i in [1..#iter] do
        iter[i](~param[i],~value[i]);

        if value[i] ne SetPartitionIterator_DONE_FLAG or i eq #iter then
          break;
        end if;

        iter[i], param[i], value[i] := Get_SetPartitionIterator_Helper(comp[i],SP,degree_bound,false);
        iter[i](~param[i],~value[i]);
      end for;

      if value[#value] eq SetPartitionIterator_DONE_FLAG then VALUE := Tilde_N_n_e_Iterator_DONE_FLAG;
      else 
        arcs := []; 
      
        for i in [1..#value] do
          cs := comp_shifts[i];
          next := SetPartition_Shift(1+cs,comp[i]+cs,cs,cs,value[i]);
          arcs := arcs cat Get_ArcSet(next);
        end for;

        VALUE := Compute_N_Lambda_e(ArcSet_to_SetPartition(arcs,1,N)) * n_Lambda_e;
        Gamma := ArcSet_to_SetPartition(arcs cat Get_ArcSet(frame),1,N);
        VALUE := VALUE * (MV ! x - 1)^(N-#Gamma);

      end if;

      PARAMS := [* iter, param, value *];

    end if;

  end procedure;

  return it, rvals, 1;  

end function;



/*
 * Function: Get_CompositionsHelper, Get_WeakCompositions, Get_Compositions
 * ------------------------------------------------------------------------
 * Get_WeakCompositions(num_parts,sum) returns the list of lists with size
 * num_parts of nonnegative integers summing to sum.
 *
 * Get_Compositions(num_parts,sum) returns the list of lists with size
 * num_parts of positive integers summing to sum.
 *
 * Get_CompositionsHelper(min_part_size,max_part_size,num_parts,sum) returns
 * the list of lists with size num_parts of nonnegative integers in the range
 * [min_part_size..max_part_size] summing to sum.
 * 
 */

forward Get_CompositionsHelper;

Get_CompositionsHelper:=function(min_part_size,max_part_size,num_parts,sum)

  if num_parts eq 0 then
    if sum eq 0 then return [[]];
    else return [];
    end if;
  end if;

  ans := [];
  for i in [min_part_size..Min(max_part_size,sum)] do
    rest := Get_CompositionsHelper(min_part_size, max_part_size, num_parts-1, sum-i);
    for w in rest do
      Append(~ans,Append(w,i));
    end for;
  end for;
  
  return ans;

end function;


Get_WeakCompositions:=function(num_parts,sum)
  return Get_CompositionsHelper(0,sum,num_parts,sum);
end function;

Get_Compositions:=function(num_parts,sum)
  return Get_CompositionsHelper(1,sum,num_parts,sum);
end function;

GetAll_Compositions:=function(sum)
  c := [];

  for k in [0..sum] do
    c := c cat Get_Compositions(k,sum);
  end for;

  return c;
end function;



/*
 * Function: Compute_Tilde_N_n_e, FormatOutput
 * ---------------------------------------------------------------------
 * This method tries to compute the function Tilde_N_{n,e}(q) defined in
 * Section 3.4 of "Combinatorial methods ...". Its four inputs are
 * n, e, SP, CrConn.  SP is the output from Get_SetPartitionUpTo(m)
 * for some m, as usual.  CrConn is the output from 
 * Get_CrConnSetPartitionsUpTo(n,SP,e); 
 *
 * The method returns three values:
 *   return_f is the polynomial in x and y, with degree e in y, such
 *            that the coefficient of y^i is Tilde_N_{n,i}(x).
 *   fns is the list of tuples <i, f> for i <= e with 
 *       f = Tilde_N_{n,i}(x) nonzero.
 *   a_vecs is the list of tuples <n,i,[a_0,a_1,...,a_ell]> corresponding
 *          to the tables in Appendix A of "Combinatorial methods ...". 
 *          These are such that Tilde_N_{n,i}(x) = Sum_k a_k (x-1)^k.  
 * 
 * If the algorithm is not successful, the outputs are 0,[],[].
 *
 * FormatOutput is a helper function for this method.
 *
 */

FormatOutput:=function(f,n,e)

  // Truncate output polynomial to be of degree e in y
  return_f := MV ! 0;
  for i in [0..e] do
    return_f := return_f + MV ! Coefficient(f,2,i) * (MV ! y)^i;
  end for;

  fns := Get_CountFns(return_f);
  a_vecs := [];

  for pair in fns do

    v := [];
    g := MV ! Evaluate(pair[2],1,MV ! x+1);
    for i in [0..MV!Degree(g)] do
      Append(~v,MV!Coefficient(g,1,i));
    end for;

    Append(~a_vecs,<n,pair[1],v>);

  end for;

  return return_f, fns, a_vecs;

end function;

  
Compute_Tilde_N_n_e:=function(n,e,SP,CrConn)

  t:=Cputime();

	
  f := MV ! 0;
  if n eq 1 then f := f + 1;
  end if;


  incr := 10;
  total := 0;
  for i in [2..n+1] do
    total := total + #CrConn[i] * Binomial(n-1,i-2);
  end for;

  printf "\nEach '.\' represents %o iterations of %o\n", incr, total;
  dot := 0;


  for i in [2..n+1] do
    for Lambda in CrConn[i] do
  
      weak := Get_WeakCompositions(i-1,n+1-i);
      for w in weak do
	
        it,PARAMS,g:=Get_Tilde_N_n_e_Iterator(Lambda,w,SP,e);
        while true do

          it(~PARAMS,~g);

          if g eq Tilde_N_n_e_Iterator_DONE_FLAG then break;
          elif g eq 0 then return 0,[],[]; 
          end if;

          f := f + g;
        end while;

        PrintProcedure(~dot,incr,total);

      end for;      

    end for;
  end for;

  print "\nTime=", Cputime(t),"\n";
  return FormatOutput(f,n,e);

end function;


/*
 * Function: Compute_Tilde_N_n_e_UpTo
 * ----------------------------------------------------------------
 * Returns the list [f_1,f_2,...,f_n] where f_n is the first output 
 * of Compute_Tilde_N_n_e(i,e,Sp,CrConn);
 *
 */

Compute_Tilde_N_n_e_UpTo:=function(n,e,SP,CrConn);

  output := [];
  printf "\n"; 

  for i in [1..n] do
    printf "Computing Tilde_N_{n,e}(q) for e = %o and n = %o . . . ", e,i;
    Append(~output, Compute_Tilde_N_n_e(i,e,SP,CrConn));
  end for;

  printf "\n";
  return output;

end function;


/*
 * Function: Compute_N_n_e_From_Tilde_N_n_e, Compute_f_e_i
 * ---------------------------------------------------------------
 * The function Compute_N_n_e_From_Tilde_N_n_e takes three inputs:
 *   n_scalar is a positive integer 
 *   e is an integer in {1,2,...,8}
 *   Tilde_N_n_e_Outputs is the (successful) output from 
 *
 *     Compute_Tilde_N_n_e_UpTo(2*e+1,e,SP,CrConn);
 *   
 *   where SP is the output from Get_SetPartitionsUpTo(m) for some m > 0
 *   and CrConn is the output from Get_CrConnSetPartitionsUpTo(2*e+1,SP). 
 *
 * The function has three outputs: 
 *   N_n_e_Specific is the polynomial is q giving  N_{n_scalar,e}(q).
 *   N_n_e_Divided is polynomial in n and q equal to N_{n,e}(q) / q^{n-e-2} 
 *     for sufficiently large n >> e.
 *   f_e_i is a list of tuples <e,i,[a_0,a_1,...,a_ell]> corresponding
 *     to the tables in Appendix B of "Combinatorial methods ...".  These are
 *     such that N_n_e_Divided = Sum_{i=1}^{2*e} c_{e,i}! / e! * (q-1)^i * f_{e,i}(n-2e-1)
 *     where f_{e,i}(x) = a_0 + a_1 x + a_2 x^2 + ...
 *
 * There is no error-checking in this method; the input Tilde_N_n_e_Outputs
 * is assumed to be correct. The outputs are computed from the inputs using 
 * Equation (3.7) and the argument preceding Theorem 3.5 in "Combinatorial
 * methods ...".
 *
 * One can use this function to check that our methods give the same 
 * functions N_{n,e}(q) for e <= 8 and n <= 13 as those listed in 
 * Isaacs's and Evseev's papers.
 *
 * Compute_f_e_i is a helper function, calculating the output f_e_i from 
 * N_n_e_Divided.
 *
 */

N_n_e_Ring<q,n>:=PolynomialRing(Rationals(),2);

Compute_f_e_i:=function(e,N_n_e_Divided)

  f_e_i := [];

  f := N_n_e_Divided * Factorial(e);
  f := N_n_e_Ring ! Evaluate(f, 1, N_n_e_Ring ! q + 1);
  f := N_n_e_Ring ! Evaluate(f, 2, N_n_e_Ring ! n + 2*e+1);

  for i in [0..N_n_e_Ring ! Degree(f,1)] do

    summand := N_n_e_Ring ! Coefficient(f,1,i) / Factorial(Round(1/2+Abs(1/2+e-i)));

    if summand ne 0 then 
      vec := [];
      for j in [0..N_n_e_Ring ! Degree(summand,2)] do
        Append(~vec, N_n_e_Ring ! Coefficient(summand,2,j));
      end for;

      Append(~f_e_i, <e,i,vec>);
    end if;
  end for;

  return f_e_i;

end function;


Compute_N_n_e_From_Tilde_N_n_e:=function(n_scalar,e,Tilde_N_n_e_Outputs)

  N_n_e_Ring<q,n>:=PolynomialRing(Rationals(),2);
  N_n_e_Divided := 0;
  N_n_e_Specific := 0;

  d_set := GetAll_Compositions(e);
  for d in d_set do
 
    for sum in [1..2*e+#d] do
  
      // The second input bounds the part size of the resulting compositions
      c_set := Get_CompositionsHelper(1, 2*Max(d)+1, #d, sum);
      for c in c_set do

        prod := 1; 
        binom := 1;

        for i in [1..#c] do

          f := 0;
          if c[i] le 2*d[i]+1 then
            f := MV ! Coefficient(Tilde_N_n_e_Outputs[c[i]],2,d[i]);
            f := N_n_e_Ring ! Evaluate(f,1, N_n_e_Ring ! q);
          end if;
          prod := prod * f;

          binom := binom * (N_n_e_Ring ! n-sum+#c-i) / (#c+1-i);
        end for;

        if prod ne 0 then 
          N_n_e_Divided := N_n_e_Divided + (N_n_e_Ring ! q)^(3*e-sum) * prod * binom;

          if sum lt n_scalar then
            N_n_e_Specific := N_n_e_Specific + N_n_e_Ring ! Evaluate((N_n_e_Ring ! q)^(n_scalar-1-sum) * prod * binom, 2, n_scalar);
          end if;
        end if;

      end for;
    end for;
  end for;

  N_n_e_Divided := N_n_e_Divided div (N_n_e_Ring ! q)^(2*e-1);
  return N_n_e_Specific, N_n_e_Divided, Compute_f_e_i(e, N_n_e_Divided);

end function;

