// brandt.magma
// andrew snowden
// 1/6/11

// use load "brandt.magma"; to get this in magma


/************
 * StructuredDecomposition
 *
 * Description:
 *   Groups the subspaces of the Brandt module according to their Hecke action.
 * Parameters:
 *   B  - A Brandt module.
 *   bd - A bound on how far to compute.
 * Return:
 *   A list, each element of which is a list of elements of Decomposition(B) which have
 *   the same Hecke action.
 */
function StructuredDecomposition(B, bd)

  // check parameters
  assert Type(B) eq ModBrdt;
  assert Type(bd) eq RngIntElt;

  // decompose module
  A:=Decomposition(B, bd);
  N:=Level(B);

  // set up loop
  bucketlist:=[A];
  p:=2;
  while p lt bd do

    // go to next prime not dividing the level
    p:=NextPrime(p);
    if N mod p eq 0 then
      continue;
    end if;

    newbucketlist:=[];
    for bucket in bucketlist do

      n:=#bucket;
      if n eq 1 then
        Append(~newbucketlist, bucket);
        continue;
      end if;

      bucket2:=bucket;
      hecke:=[CharacteristicPolynomial(HeckeOperator(V,p)) : V in bucket];
      ParallelSort(~hecke, ~bucket2);
      i:=1;
      while i le n do
        P:=hecke[i];
        newbucket:=[];
        while (i le n) and (hecke[i] eq P) do
          Append(~newbucket, bucket2[i]);
          i:=i+1;
        end while;
        Append(~newbucketlist, newbucket);
      end while;

    end for;

    bucketlist:=newbucketlist;

  end while;

  return bucketlist;

end function;



/************
 * NewDecomposition
 *
 * Description:
 *   Computes the list of new spaces of a Brandt module.
 * Parameters:
 *   B  - A Brandt module.
 *   bd - A bound on how far to compute.
 * Return:
 *   The list of elements of Decomposition(B) which are new.
 */
function NewDecomposition(B, bd)

  // decompose module
  A:=StructuredDecomposition(B, bd);

  // find new ones
  list:=[];
  for bucket in A do
    if #bucket eq 1 then
      Append(~list, bucket[1]);
    end if;
  end for;

  // return what we found
  return list;

end function;



/************
 * MultiStructuredDecomposition
 *
 * Description:
 *   Groups the subspaces of multiple Brandt modules according to their Hecke action.
 * Parameters:
 *   B - A sequence of Brandt modules.
 *   N - Bad primes (don't consider T_p with p dividing N).  Must include all level.
 *   bd - A bound on how far to compute.
 * Return:
 *   A sequence, each item of which is a sequence of the form [x_1, ..., x_n], where
 *   n is the number of Brandt modules in B.  Each x_i is a sequence of elements
 *   of Decomposition(B[i]) which have the same Hecke action.  Furthermore, the
 *   elements of x_i and of x_j also have the same Hecke action.
 */
function MultiStructuredDecomposition(B, N, bd)

  // check parameters
  n:=#B;
  for i in [1..n] do
    assert Type(B[i]) eq ModBrdt;
    assert N mod Level(B[i]) eq 0;
  end for;
  assert Type(N) eq RngIntElt;
  assert Type(bd) eq RngIntElt;

  // build a big table with everything we need
  table:=[];
  primes:=[p : p in [1..bd] | IsPrime(p) and (N mod p ne 0)];
  for i in [1..n] do
    A:=Decomposition(B[i], bd);
    for U in A do
      h:=[CharacteristicPolynomial(HeckeOperator(U,p)) : p in primes];
      Append(~table, [* i, U, h *]);
    end for;
  end for;

  // start grouping
  ret:=[];
  count:=#table;
  while count ne 0 do

    // get first hecke list
    h:=table[1][3];

    // now find everything with hecke entry h
    list:=[[] : i in [1..n]];
    i:=1;
    while i le count do
      if table[i][3] ne h then
        i:=i+1;
        continue;
      end if;
      j:=table[i][1];
      Append(~list[j], table[i][2]);
      Remove(~table, i);
      count:=count-1;
    end while;
    Append(~ret, list);

  end while;

  // return
  return ret;

end function;



/************
 * MultiNewDecomposition
 *
 * Description:
 *   Like MultiStructuredDecomposition, but only returns new spaces.
 * Parameters:
 *   B - A sequence of Brandt modules.
 *   N - Bad primes (don't consider T_p with p dividing N).  Must include all level.
 *   bd - A bound on how far to compute.
 * Return:
 *   A sequence, each item of which is a sequence of the form [x_1, ..., x_n], where
 *   n is the number of Brandt modules in Blist.  Each x_i is a new subspace of B[i]
 *   with the same Hecke action.
 */
function MultiNewDecomposition(B, N, bd)
  A:=MultiStructuredDecomposition(B, N, bd);
  ret:=[];
  for list in A do
    list2:=[];
    flag:=true;
    for x in list do
      if #x ne 1 then
        flag:=false;
        break;
      end if;
      Append(~list2, x[1]);
    end for;
    if flag then
      Append(~ret, list2);
    end if;
  end for;
  return ret;
end function;



/************
 * FindSubspaceGeneral
 *
 * Description:
 *   Finds a subspace of a Brandt module with specified Hecke action.
 * Parameters:
 *   B - A Brandt module
 *   L - A list of pairs (p, P) where p is a prime and P is a polynomial
 *       of degree d over the integers
 * Return:
 *   The list of elements of Decomposition(B) for which the Hecke operator T_p has
 *   characteristic polynomial P for those p in the list L.
 */
function FindSubspaceGeneral(B, L)

  // check parameters
  assert Type(B) eq ModBrdt;
  assert Type(L) eq SeqEnum;
  first:=true;
  for pair in L do
    assert Type(pair) eq List;
    assert #pair eq 2;
    p:=pair[1];
    P:=pair[2];
    assert Type(p) eq RngIntElt;
    assert IsPrime(p);
    assert Type(P) eq RngUPolElt;
    if first then
      d:=Degree(P);
      first:=false;
    else
      assert Degree(P) eq d;
    end if;
  end for;

  // decompose the brandt module
  A:=Decomposition(B, 100);

  // examine each space in turn
  ret:=[];
  for U in A do
    if Dimension(U) eq d then
      flag:=true;
      for pair in L do
        p:=pair[1];
        P:=pair[2];
        T:=HeckeOperator(U, p);
        if CharacteristicPolynomial(T) ne P then
          flag:=false;
          break;
        end if;
      end for;
      if flag then
        Append(~ret, U);
      end if;
    end if;
  end for;

  // return what we found
  return ret;

end function;



/************
 * FindSubspaceRational
 *
 * Description:
 *   Same as FindSubspaceGeneral but only applies when the coefficient field is Q.
 *   The parameters are simpler in this case.
 * Parameters:
 *   B - A Brandt module
 *   L - A list of pairs (p, e) where p is a prime and e is an integer specifying the
 *       eigenvalue of T_p
 * Return:
 *   The list of elements of Decomposition(B) d for which the Hecke operator T_p has
 *   characteristic polynomial P for those p in the list L.
 */
function FindSubspaceRational(B,L)

  // check parameters
  assert Type(B) eq ModBrdt;
  assert Type(L) eq SeqEnum;
  for pair in L do
    assert Type(pair) eq List;
    assert #pair eq 2;
    p:=pair[1];
    e:=pair[2];
    assert Type(p) eq RngIntElt;
    assert IsPrime(p);
    assert Type(e) eq RngIntElt;
  end for;

  // decompose the brandt module
  A:=Decomposition(B, 100);

  // examine each space in turn
  ret:=[];
  for U in A do
    if Dimension(U) eq 1 then
      flag:=true;
      for pair in L do
        p:=pair[1];
        e:=pair[2];
        T:=HeckeOperator(U, p);
        if T[0][0] ne e then
          flag:=false;
          break;
        end if;
      end for;
      if flag then
        Append(~ret, U);
      end if;
    end if;
  end for;

  // return what we found
  return ret;

end function;



/************
 * ComputeNormGeneral
 *
 * Description:
 *   Computes the norm of a specified newform in a Brandt module.  The newform
 *   is specified as the eigenvector of T_p with eigenvalues e in the space V.
 * Parameters:
 *   V - An irreducible piece of a Brandt module.
 *   K - The coefficient field.
 *   p - A prime number.
 *   e - An element of K.
 * Return:
 *   The norm of the specified newform, given as an ideal of the maximal order of K.
 * Notes:
 *   V can be found by using FindSubspaceGeneral.
 */
function ComputeNormGeneral(V, K, p, e)

  // compute a basis for V
  assert Type(V) eq ModBrdt;
  basis:=Basis(V);

  // handle the rational case
  d:=Degree(K);
  if d eq 1 then
    v:=basis[1];
    return Norm(v);
  end if;

  // do a bit more error checking
  assert Type(K) eq FldNum;
  assert Parent(e) eq K;
  assert Type(p) eq RngIntElt;
  assert IsPrime(p);

  // find eigenvector
  T:=HeckeOperator(V, p);
  TK:=ChangeRing(T, K);
  E:=Eigenspace(TK, e);
  assert Dimension(E) eq 1;
  v:=Basis(E)[1];

  // compute the norm of v
  //   this is slightly annoying:  V is an inner product space over the rationals, and
  //     we'd like to regard v as an element of V \otimes K.  i don't know how to do this
  //     well in magma.
  //   instead, we create the matrix X of inner products of the basis "basis" and then
  //     create an inner product space over K with this basis.
  //   we then create an element of this new space from our old element v and take its norm
  X:=Matrix(K,d,d,[[InnerProduct(basis[i], basis[j]) : i in [1..d]] : j in [1..d]]);
  V:=VectorSpace(K, d, X);
  f:=V![v[i] : i in [1..d]];
  norm:=Norm(f);

  // scale the norm to account for v not being normalized
  //   the basis "basis" is an integral basis
  //   v is the eigenvector, we treat it as a tuple of coefficients in the basis "basis"
  //   v is often not integral, so we must scale it
  //   to do this, we look at the fractional ideal generated by the coefficient of v,
  //     find a generator, and scale v by this.  (this only works if that ideal is principal.)
  O:=MaximalOrder(K);
  I:=ideal<O|[v[i] : i in [1..d]]>;
  return I^(-2)*norm;

end function;



/************
 * ComputeNormRational
 *
 * Description:
 *   Same as ComputeNormGeneral but only applies when the coefficient field is Q.
 *   The parameters are simpler in this case.
 * Parameters:
 *   V - An irreducible piece of a Brandt module.
 * Return:
 *   The norm of the specified newform, an integer.
 * Notes:
 *   V can be found by using FindSubspaceRational.
 */
function ComputeNormRational(V)
  assert Type(V) eq ModBrdt;
  v:=Basis(V)[1];
  return Norm(v);
end function;



/************
 * ComputeNorms
 *
 * Description:
 *   Given several corresponding Brandt modules, simulateneously compute their norms.
 * Parameters:
 *   V - A list of irreducible subspaces of Brandt modules which correspond.
 * Return:
 *   A list [* K, p, e, x, y *] where K is the coefficient field, p is a prime, e is an
 *   eigenvalue of T_p which generates K, x_i is the norm of the normalized element of V_i
 *   on which T_p acts by e, and y_i is the norm from K to Q of x_i.
 */
function ComputeNorms(V)

  // check parameter
  n:=#V;
  for i in [1..n] do
    assert Type(V[i]) eq ModBrdt;
  end for;
  N:=Level(V[1]);
  d:=Dimension(V[1]);
  for i in [1..n] do
    assert Level(V[i]) eq N;
    assert Dimension(V[i]) eq d;
  end for;

  // find a prime not dividing the level such that the char poly of T_p is irreducible
  p:=3;
  found:=false;
  while not found do
    assert p le 5000;
    p:=NextPrime(p);
    if N mod p eq 0 then
      continue;
    end if;
    T:=HeckeOperator(V[1], p);
    P:=CharacteristicPolynomial(T);
    found:=IsIrreducible(P);
  end while;

  // build number field and pick eigenvalue of T_p
  K:=NumberField(P);
  TK:=ChangeRing(T, K);
  if Evaluate(CharacteristicPolynomial(TK), K.1) eq 0 then
    e:=K.1;
  else
    e:=Random(Eigenvalues(TK))[1];
  end if;

  // compute norms
  x:=[];
  y:=[];
  for i in [1..n] do
    x[i]:=ComputeNormGeneral(V[i], K, p, e);
    y[i]:=IntegerRing()!Norm(x[i]);
  end for;

  // return
  return [* K, p, e, x, y *];

end function;



/************
 * DisplayNorms
 *
 * Description:
 *   Compute and print norms of corresponding Brandt modules.
 * Parameters:
 *   fp - A file to write to.
 *   V - A list of irreducible subspaces of Brandt modules which correspond.
 */
procedure DisplayNorms(fp, V)

  // compute a bunch of stuff
  list:=ComputeNorms(V);
  K:=list[1];
  p:=list[2];
  e:=list[3];
  x:=list[4];
  y:=list[5];
  n:=#V;
  N:=Level(V[1]);
  d:=Dimension(V[1]);
  disc:=[PrimeDivisors(Discriminant(V[i])) : i in [1..n]];
  cond:=[PrimeDivisors(Conductor(V[i])) : i in [1..n]];
  O:=MaximalOrder(K);
  assert Degree(K) eq d;

  // a break line
  fprintf fp, "\n=========================================\n\n";

  // basic information about the spaces we're looking at
  if d eq 1 then
    fprintf fp, "Coefficient field is the rationals.\n";
    fprintf fp, "Eigenvalues of small Hecke operators:\n";
    for p in [2..31] do
      if IsPrime(p) and (N mod p ne 0) then
        T:=HeckeOperator(V[1], p);
        fprintf fp, "T_%o : %o\n", p, T[1][1];
      end if;
    end for;
  else
    fprintf fp, "Coefficient field: %o\n", K;
    fprintf fp, "Generated by T_%o, eigenvalue %o\n", p, e;
    fprintf fp, "Basis of maximal order: ";
    b:=Basis(O);
    assert #b eq d;
    for i in [1..d] do
      fprintf fp, "%o", K!b[i];
      if i ne d then
        fprintf fp, ", ";
      end if;
    end for;
    fprintf fp, "\n";
    fprintf fp, "Characteristic polynomials of small Hecke operators:\n";
    for p in [2..31] do
      if IsPrime(p) and (N mod p ne 0) then
        P:=CharacteristicPolynomial(HeckeOperator(V[1], p));
        fprintf fp, "T_%o : %o\n", p, P;
      end if;
    end for;
  end if;
  fprintf fp, "\n";

  // a subprocedure to display factorizations of ideals
  procedure DisplayFactorization(I)
    S:=Factorization(I);
    for z in S do
      P:=z[1];
      e:=z[2];
      p:=Norm(P);
      g1:=Basis(P)[1];
      g2:=Basis(P)[2];
      fprintf fp, "   < p=%o, e=%o, g1=%o, g2=%o >\n", p, e, g1, g2;
    end for;
  end procedure;

  // show the norms
  if d gt 1 then
    fprintf fp, "Norms, as ideals of the coefficient field:\n";
    for i in [1..n] do
      fprintf fp, "%o %o :\n", disc[i], cond[i];
      DisplayFactorization(x[i]);
    end for;
    fprintf fp, "\n";
  end if;

  // show the norms
  if d gt 1 then
    fprintf fp, "Norms down to Q:\n";
  else
    fprintf fp, "Norms:\n";
  end if;
  for i in [1..n] do
    fprintf fp, "%o %o : %o\n", disc[i], cond[i], Factorization(y[i]);
  end for;
  fprintf fp, "\n";

  // flush
  Flush(fp);

end procedure;


function ComputeBrandtModules(primes)

/*
  p1:=primes[1];
  p2:=primes[2];
  p3:=primes[3];
  p4:=primes[4];

  print "Computing Brandt module 1/8";
  B1:=BrandtModule(p1, p2*p3*p4);
  print "Computing Brandt module 2/8";
  B2:=BrandtModule(p2*p3*p4, p1);
  print "Computing Brandt module 3/8";
  B3:=BrandtModule(p2, p1*p3*p4);
  print "Computing Brandt module 4/8";
  B4:=BrandtModule(p1*p3*p4, p2);
  print "Computing Brandt module 5/8";
  B5:=BrandtModule(p3, p1*p2*p4);
  print "Computing Brandt module 6/8";
  B6:=BrandtModule(p1*p2*p4, p3);
  print "Computing Brandt module 7/8";
  B7:=BrandtModule(p4, p1*p2*p3);
  print "Computing Brandt module 8/8";
  B8:=BrandtModule(p1*p2*p3, p4);

  return [B1,B2,B3,B4,B5,B6,B7,B8];

*/


  // get primes
  p1:=primes[1];
  p2:=primes[2];
  p3:=primes[3];
  p4:=primes[4];
  p5:=primes[5];

  // compute brandt modules
  print "Computing Brandt module 1/16";
  B1:=BrandtModule(p1, p2*p3*p4*p5);
  print "Computing Brandt module 2/16";
  B2:=BrandtModule(p2, p1*p3*p4*p5);
  print "Computing Brandt module 3/16";
  B3:=BrandtModule(p3, p1*p2*p4*p5);
  print "Computing Brandt module 4/16";
  B4:=BrandtModule(p4, p1*p2*p3*p5);
  print "Computing Brandt module 5/16";
  B5:=BrandtModule(p5, p1*p2*p3*p4);
  print "Computing Brandt module 6/16";
  B6:=BrandtModule(p1*p2*p3, p4*p5);
  print "Computing Brandt module 7/16";
  B7:=BrandtModule(p1*p2*p4, p3*p5);
  print "Computing Brandt module 8/16";
  B8:=BrandtModule(p1*p2*p5, p3*p4);
  print "Computing Brandt module 9/16";
  B9:=BrandtModule(p1*p3*p4, p2*p5);
  print "Computing Brandt module 10/16";
  B10:=BrandtModule(p1*p3*p5, p2*p4);
  print "Computing Brandt module 11/16";
  B11:=BrandtModule(p1*p4*p5, p2*p3);
  print "Computing Brandt module 12/16";
  B12:=BrandtModule(p2*p3*p4, p1*p5);
  print "Computing Brandt module 13/16";
  B13:=BrandtModule(p2*p3*p5, p1*p4);
  print "Computing Brandt module 14/16";
  B14:=BrandtModule(p2*p4*p5, p1*p3);
  print "Computing Brandt module 15/16";
  B15:=BrandtModule(p3*p4*p5, p1*p2);
  print "Computing Brandt module 16/16";
  B16:=BrandtModule(p1*p2*p3*p4*p5, 1);

  // make list of brandt modules
  return [B1,B2,B3,B4,B5,B6,B7,B8,B9,B10,B11,B12,B13,B14,B15,B16];

end function;



/************
 * DoEverything
 *
 * Description:
 *   Compute Brandt modules, find corresponding new subspaces, compute norms and output to file.
 * Parameters:
 *   file - Filename for output
 *   primes - a list of four odd primes
 */
procedure DoEverything(file, B)

  // compute brandt modules
  //B:=ComputeBrandtModules(primes);

  // decompose
  print "Decomposing Brandt modules";
  A:=MultiNewDecomposition(B, 3*5*7*11*13, 35);

  // output
  print "Computing norms";
  fp:=Open(file, "wt");
  for V in A do
    DisplayNorms(fp, V);
  end for;
  delete fp;

  print "Done!";

end procedure;