// This Magma script is for computing models of modular curves
// in a relative way. This is version 10. It is a slight modification
// of version 9 where a few steps are added in the functions
// ratfuncrep and genusoneratfuncrep to ensure that the answer obtained
// is unique. 

if (not assigned m) then
  printf "This script requires that m be passed as a command line argument.\n";
  printf "The usage is something like magma m:=5 models9.txt";
  quit;
else
  subnum := StringToInteger(m);
end if;

// These lines are to fool the Magma parser into working.
// These are identifiers that are used inside if statements
// that are referred to again later inside other if statements

E2 := 0;
xmodfunc := 0;
newC := 0;
ymodfunc := 0;
AA := 0;
BB := 0;
C2 := 0;
s2 := 0;
phi3 := 0;
newratfunc := 0;
haupt := 0;
newhaupt := 0;

// Step 1 - Build covering group and read the stored Fourier expansion.

load "gl2data.txt";
// Select a subgroup
logfilename := "model" cat IntegerToString(subnum) cat ".out";
System("rm " cat logfilename);
SetLogFile(logfilename);
N := newsublist[subnum][2];
sub := newsublist[subnum][3];
cover := newsublist[subnum][7];
printf "Working on X_%o. It covers X_%o.\n",subnum,cover;
ind := newsublist[subnum][1] div newsublist[cover][1];
printf "The degree of that covering map is %o.\n",ind;
covergpset := { newsublist[i][7] : i in [1..#newsublist] };
K<zeta> := CyclotomicField(N);
R<q> := PuiseuxSeriesRing(K);
// We're going with this for now, and we'll see if we need to modify it
// down the road. 
// prec := 110*N;
prec := 110*N;

if subnum eq 1 then
  gencover := 0;
  else
  gencover := newsublist[cover][5];
end if;

// Read the data from the stored files

// This function returns true if a file exists, and false if not
function fileexists(str)
  retval := System("test -e " cat str);
  retval := retval div 256;
  ret := false;
  if (retval eq 0) then
    ret := true;
  end if;
  return ret;  
end function;

dontdoset := {};
computedset := {};
if not fileexists("doneset.txt") then
  PrintFile("doneset.txt","dontdoset := {};");
  PrintFile("doneset.txt","computedset := {};");
  PrintFile("doneset.txt","return dontdoset, computedset;");
else
  str := Read("doneset.txt");
  dontdoset, computedset := eval str;
  if (subnum in dontdoset) then
    printf "This curve covers one with only finitely many rational points.\n";
    printf "We don't need to worry about it!\n";
    quit; 
  end if;
end if;

// If gencover eq 0 then we have stored the Fourier expansion of a hauptmodul
// We suppose that the Fourier expansion was stored with z as the root of unity.
// This root of unity may be different than zeta. The power series ring
// should be defined in the file. The modular function should be called Haupt

if gencover eq 0 then
  filename := "modfunc" cat IntegerToString(cover) cat ".txt";
  str := Read(filename);
  haupt := eval str;
  cycfield := BaseRing(Parent(haupt));
  ordofz := 2*Degree(cycfield);
  pow := Floor(N/ordofz);
  expdenom := ExponentDenominator(haupt);
  newhaupt := R!0;
  for n in [Valuation(haupt)*expdenom..AbsolutePrecision(haupt)*expdenom-1] do
    cof := Eltseq(Coefficient(haupt,n/expdenom));
    term := R!(&+[ cof[i]*zeta^(pow*(i-1)) : i in [1..#cof]]*q^(n/expdenom));
    newhaupt := newhaupt + term;
  end for;
  newhaupt := newhaupt + BigO(q^AbsolutePrecision(haupt));
  haupt := newhaupt;
end if;

// If gencover eq 1 then we have stored the Fourier expansions of the
// x and y coordinate functions on the elliptic curve, as well as its
// model.

xcoord := 0;
ycoord := 0;
expdenomx := 0;
expdenomy := 0;
alist := [0,0,0,0,0];
if gencover eq 1 then
  filename := "eq" cat IntegerToString(cover) cat ".txt";
  str := Read(filename);
  ee := eval str;
  alist := aInvariants(ee);
  filename := "modfunc" cat IntegerToString(cover) cat ".txt";
  str := Read(filename);
  xcoord, ycoord := eval str;
  cycfield := BaseRing(Parent(xcoord));
  ordofz := 2*Degree(cycfield);
  pow := Floor(N/ordofz);
  newxcoord := R!BigO(q^AbsolutePrecision(xcoord));
  newycoord := R!BigO(q^AbsolutePrecision(ycoord));
  expdenomx := ExponentDenominator(xcoord);
  expdenomy := ExponentDenominator(ycoord);
  for n in [Valuation(xcoord)*expdenomx..AbsolutePrecision(xcoord)*expdenomx-1] do
    cof := Eltseq(Coefficient(xcoord,n/expdenomx));
    term := R!(&+[ cof[i]*zeta^(pow*(i-1)) : i in [1..#cof]]*q^(n/expdenomx));
    newxcoord := newxcoord + term;
  end for;
  for n in [Valuation(ycoord)*expdenomy..AbsolutePrecision(ycoord)*expdenomy-1] do
    cof := Eltseq(Coefficient(ycoord,n/expdenomy));
    term := R!(&+[ cof[i]*zeta^(pow*(i-1)) : i in [1..#cof]]*q^(n/expdenomy));
    newycoord := newycoord + term;
  end for;
  xcoord := newxcoord;
  ycoord := newycoord;
end if;

// This function takes a subgroup G of GL(2,Z/2^nZ) and lifts
// it to a subgroup of GL(2,Z/2^mZ) for m > n.

function liftsub(G,n,m)
  H2 := GL(2,Integers(2^m));
  H3 := GL(2,Integers(2^n));
  genlist := [];
  for g in Generators(G) do
    Append(~genlist,H2!g);
  end for;
  phi := hom<H2->H3 | [H3!H2.i : i in [1..#Generators(H2)]]>;
  ker := Kernel(phi);
  // Add generators for the kernel of GL_2(Z/2^m Z) -> GL_2(Z/2^n Z)
  return sub<H2 | genlist, ker>;
end function;

// This function takes a subgroup G of GL(2,Integers(N)) and
// compute the genus of the corresponding modular curve, as
// well as the index, number of cusps, e2 and e3.

function genus2(G)
  md := Modulus(BaseRing(G));
  H := SL(2,Integers(md));
  S := H![0,-1,1,0];
  T := H![1,1,0,1];
  phi, perm := CosetAction(H,G meet H);
  lst := [phi(S),phi(T),phi(S*T)];
  //printf "Permutation for S = %o.\n",phi(S);
  //printf "Permutation for T = %o.\n",phi(T);
  //printf "Permutation for S*T = %o.\n",phi(S*T);
  cs := [CycleStructure(lst[i]) : i in [1..3]];
  gen := -2*Degree(perm) + 2;
  einfty := #Orbits(sub<perm | lst[2]>);
  e2 := #Fix(lst[1]);
  e3 := #Fix(lst[3]);
  ind := Degree(perm);
  for j in [1..3] do
    for i in [1..#cs[j]] do
      gen := gen + (cs[j][i][1]-1)*cs[j][i][2];
    end for;
  end for;
  gen := gen div 2;
  //printf "The genus = %o.\n",gen;
  genhur := 1 + (ind/12) - (e2/4) - (e3/3) - (einfty/2);
  //printf "The Hurwitz formula is %o = 1 + %o/12 - %o/4 - %o/3 - %o/2.\n",genhur,ind,e2,e3,einfty;
  return gen, ind, einfty, e2, e3;
end function;


// This function takes an element r of (Z/NZ) and
// lifts it to an integer between 1 and N

function lift(r,N)
  liftr := Integers()!r;
  liftr := liftr mod N;
  if (liftr eq 0) then
    liftr := N;
  end if;
  return liftr;
end function;

// Simplification function. This takes as input a rational function f
// and returns a simplier version of f, obtained by composing
// with a linear fractional transformation. This transformation is also
// returned.

/*
function size(f)
  num := Numerator(f);
  denom := Denominator(f);
  coflist := [ Denominator(Coefficient(num,n)) : n in [0..Degree(num)]];
  coflist := coflist cat [ Denominator(Coefficient(denom,n)) : n in [0..Degree(denom)]];
  mult := LCM(coflist);
  num := num*mult;
  denom := denom*mult;
  if Degree(num) gt 0 then
    ret := &+[ Log(1+AbsoluteValue(Numerator(Coefficient(num,n)))) : n in [0..Degree(num)]];
  else
    ret := Log(1);
  end if;
  if Degree(denom) gt 0 then
    ret := ret + (&+[ Log(1+AbsoluteValue(Numerator(Coefficient(denom,n)))) : n in [0..Degree(denom)]]);
  end if;
  return ret;
end function;
*/

function size(f)
  return #Sprintf("%o",f);
end function;

// Supersimplify functions - This is probably only necessary on
// functions coming from index 3 or 4 covers because the rational
// functions that are output are ridiculously complicated

function transcheck(f)
  //printf "Finding optimal translation.\n";
  t := Parent(f).1;
  done := true;
  if size(Evaluate(f,t-1)) lt size(f) then
    upp := 0;
    mid := -1;
    low := -2;
    while size(Evaluate(f,t+low)) lt size(Evaluate(f,t+mid)) do
      upp := mid;
      mid := low;
      low := low*2;
    end while;
    done := false;
  end if;
  if size(Evaluate(f,t+1)) lt size(f) then
    low := 0;
    mid := 1;
    upp := 2;
    while size(Evaluate(f,t+2*mid)) lt size(f) do
      low := mid;
      mid := 2*mid;
      upp := 2*mid;
    end while;
    done := false;
  end if;
  if done eq false then
    uppsiz := size(Evaluate(f,t+upp));
    midsiz := size(Evaluate(f,t+mid));
    lowsiz := size(Evaluate(f,t+low));
    round := 0;
    while done eq false do
      round := round + 1;
      //printf "Bisection method, round %o.\n",round;
      //printf "low = %o, lowsiz = %o.\n",low,lowsiz;
      //printf "mid = %o, midsiz = %o.\n",mid,midsiz;
      //printf "upp = %o, uppsiz = %o.\n",upp,uppsiz;
      if (upp-low) le 1 then
        done := true;
      else
        check1 := Round((low+mid)/2);
        check2 := Round((mid+upp)/2);
        newsiz1 := size(Evaluate(f,t+check1));
        newsiz2 := size(Evaluate(f,t+check2));
        sizelist := [lowsiz,newsiz1,midsiz,newsiz2,uppsiz];
        ind := Index(sizelist,Min(sizelist));
        if (ind eq 2) then
          upp := mid;
          uppsiz := midsiz;
          mid := check1;
          midsiz := newsiz1;
        end if;
        if (ind eq 3) then
          low := check1;
          lowsiz := newsiz1;
          upp := check2;
          uppsiz := newsiz2;
        end if;
        if (ind eq 4) then
          low := mid;
          lowsiz := midsiz;
          mid := check2;
          midsiz := newsiz2;
        end if;
      end if;
    end while;
    f2 := Evaluate(f,t+mid);
    if size(f2) lt size(f) then
      //printf "It is %o.\n",mid;
      retM := Matrix([[1,mid],[0,1]]);
      //printf "It is %o.\n",mid;
      retM := Matrix([[1,mid],[0,1]]);
    else
      //printf "It is %o.\n",0;
      f2 := f;
      retM := Matrix([[1,0],[0,1]]);
    end if;
  else
    f2 := f;
    //printf "It is %o.\n",0;
    retM := Matrix([[1,0],[0,1]]);
  end if;
  return f2,retM;
end function;

// This function returns the optimal scaling of the
// polynomial g with relatively prime integer coefficients.

function scale(g)
  coflist := Coefficients(g);
  ret := 1;
  if #[ i : i in [1..#coflist] | coflist[i] ne 0] ge 2 then
    gcd1 := GCD([Coefficient(g,i) : i in [1..Degree(g)]]);
    //printf "GCD1 = %o.\n",gcd1;
    gcd2 := GCD([Coefficient(g,i) : i in [0..Degree(g)-1]]);
    //printf "GCD2 = %o.\n",gcd2;
    //printf "Computing prime factorization.\n";
    primelist := PrimeDivisors(LCM(gcd1,gcd2));
    //printf "Done!\n";
    for j in [1..#primelist] do
      p := primelist[j];
      vallist := [ Valuation(Coefficient(g,n),p) : n in [0..Degree(g)]];
      //printf "Checking p = %o.\n",p;
      if Valuation(Coefficient(g,Degree(g)),p) gt 0 then
        //printf "Exponent should be negative.\n";
        //printf "List = %o.\n",[ vallist[n+1]/n : n in [1..Degree(g)]];
        ex := Floor(Min([ vallist[n+1]/n : n in [1..Degree(g)]]));
        ex := -ex;
      else
        //printf "Exponent should be positive.\n";
        //printf "List = %o.\n",[ vallist[n+1]/(Degree(g)-n) : n in [0..Degree(g)-1]];
        ex := Floor(Min([ vallist[n+1]/(Degree(g)-n) : n in [0..Degree(g)-1]]));
      end if;

      ret := ret*p^(ex);
    end for;
  end if;
  return ret;
end function;

function scale2(f)
  t := Parent(f).1;
  num := Numerator(f);
  denom := Denominator(f);
  coflist := [ Denominator(Coefficient(num,n)) : n in [0..Degree(num)]];
  coflist := coflist cat [ Denominator(Coefficient(denom,n)) : n in [0..Degree(denom)]];
  mult := LCM(coflist);
  num := PolynomialRing(Integers())!(num*mult);
  denom := PolynomialRing(Integers())!(denom*mult);
  //printf "Current size = %o.\n",size(f);
  //printf "Scaling numerator.\n";
  r1 := scale(num div Content(num));
  //printf "Scaling denominator.\n";
  r2 := scale(denom div Content(denom));
  primelist := PrimeDivisors(Numerator(r1)*Denominator(r1)*Numerator(r2)*Denominator(r2));
  primevals := [];
  for i in [1..#primelist] do
    v1 := Valuation(r1,primelist[i]);
    v2 := Valuation(r2,primelist[i]);
    if (v1 lt 0) and (v2 lt 0) then
      Append(~primevals,Max(v1,v2));
    else
      if (v1 gt 0) and (v2 gt 0) then
        Append(~primevals,Min(v1,v2));
      else
        Append(~primevals,0);
      end if;
    end if;
  end for;
  if #primelist gt 0 then
    scalfac := &*[ primelist[i]^primevals[i] : i in [1..#primelist]];
  else
    scalfac := 1;
  end if;
  //printf "Scaling factor %o.\n",scalfac;
  newf := Evaluate(f,scalfac*t);
  retsize := size(newf);
  //printf "Size of scaled f = %o.\n",retsize;
  return scalfac,retsize;
end function;

// Supersimplify function

function supersimplify(f)
  printf "Call to supersimplify.\n";
  t := Parent(f).1;
  j := f;
  changemat := Matrix([[1,0],[0,1]]);
  prevsize := size(j);
  prevj := j;
  prevmat := changemat;
  alldone := false;
  while alldone eq false do
  printf "Entering simplification iteration. Current size = %o.\n",prevsize;
  scal, newsize := scale2(j);
  jnew := Evaluate(j,t*scal);
  changemat := changemat*Matrix([[scal,0],[0,1]]);
  printf "Size after scaling = %o.\n",newsize;

  // Do translations - do this by
  // factoring num and denom modulo prime
  // powers

  //printf "Translation check.\n";
  num := PolynomialRing(Rationals())!Numerator(jnew);
  denom := PolynomialRing(Rationals())!Denominator(jnew);
  coflist := [ Denominator(Coefficient(num,n)) : n in [0..Degree(num)]];
  num := num*LCM(coflist);
  num := PolynomialRing(Integers())!num;
  coflist := [ Denominator(Coefficient(denom,n)) : n in [0..Degree(denom)]];
  denom := denom*LCM(coflist);
  denom := PolynomialRing(Integers())!denom;
  if Degree(num) gt 0 then
    FF := Factorization(num);
    sqrfreenum := &*[ FF[i][1] : i in [1..#FF]];
  else
    sqrfreenum := PolynomialRing(Integers())!1;
  end if;
  if Degree(denom) gt 0 then
    FF := Factorization(denom);
    sqrfreedenom := &*[ FF[i][1] : i in [1..#FF]];
  else
    sqrfreedenom := PolynomialRing(Integers())!1;
  end if;
  if (Degree(sqrfreenum) gt 0) and (Degree(sqrfreedenom) gt 0) then
    D := GCD(Discriminant(sqrfreenum),Discriminant(sqrfreedenom));
  end if;
  if Degree(sqrfreenum) eq 0 then
    D := Discriminant(sqrfreedenom);
  end if;
  if Degree(sqrfreedenom) eq 0 then
    D := Discriminant(sqrfreenum);
  end if;
  plist := PrimeDivisors(D);
  for n in [1..#plist] do
    p := plist[n];
    F := GF(p);
    RRRR<zzzz> := PolynomialRing(F);
    done := false;
    while done eq false do
      //printf "Translation check at p = %o.\n",p;
      num := PolynomialRing(Rationals())!Numerator(jnew);
      denom := PolynomialRing(Rationals())!Denominator(jnew);
      coflist := [ Denominator(Coefficient(num,n)) : n in [0..Degree(num)]];
      num := num*LCM(coflist);
      num := RRRR!num;
      coflist := [ Denominator(Coefficient(denom,n)) : n in [0..Degree(denom)]];
      denom := denom*LCM(coflist);
      denom := RRRR!denom;
      prod := RRRR!1;
      if (num ne 0) then
        prod := prod*num;
      end if;
      if (denom ne 0) then
        prod := prod*denom;
      end if;
      fac := Factorization(prod);
      //printf "Factorization = %o.\n",fac;
      if (#fac eq 1) and Degree(fac[1][1]) eq 1 then
        r := Integers()!(Roots(prod)[1][1]);
        chk1 := Evaluate(jnew,p*t+r);
        chk2 := Evaluate(jnew,p*t+(-p+r));
        size1 := size(chk1);
        size2 := size(chk2);
       //printf "Possible new sizes are %o and %o.\n",size1,size2;
        minsize := Min(size1,size2);
        if (minsize lt newsize) then
          if (size1 eq minsize) then
            jnew := chk1;
            changemat := changemat*Matrix([[p,r],[0,1]]);
            newsize := size1;
            //printf "Transformation %o. New size = %o.\n",p*t+r,newsize;
          else
            jnew := chk2;
            changemat := changemat*Matrix([[p,r-p],[0,1]]);
            newsize := size2;
            //printf "Transformation %o. New size = %o.\n",p*t+(-p+r),newsize;
          end if;
        else
          done := true;
        end if;
      else
        done := true;
      end if;
    end while;
  end for;
  printf "Translation check done. New size = %o.\n",newsize;
  j := jnew;

  // Do some rounds of random substitutions
  done := false;
  bound := 5;
  it := 0;
  ptlist := [];
  for aa in [-bound..bound] do
    for bb in [-bound..bound] do
      for cc in [-bound..bound] do
        for dd in [0..bound] do
          if GCD([aa,bb,cc,dd]) eq 1 then
            if aa*dd - bb*cc ne 0 then
              Append(~ptlist,<aa,bb,cc,dd>);
            end if;
          end if;
        end for;
      end for;
    end for;
  end for;
  printf "Doing up to 5 rounds of substitutions.\n";
  while done eq false do
    it := it + 1;
    cursize := size(j);
    //printf "Beginning iteration %o. Size = %o.\n",it,cursize;
    minsize := cursize;
    for n in [1..#ptlist] do
      pt := ptlist[n];
      aa := pt[1];
      bb := pt[2];
      cc := pt[3];
      dd := pt[4];
      jnew := Evaluate(j,(aa*t+bb)/(cc*t+dd));
      chksize := size(jnew);
      if (chksize lt minsize) then
        //printf "Index %o has size %o. pt = %o\n",n,chksize,pt;
        minsize := chksize;
        ind := n;
      end if;
    end for;
    if minsize lt cursize then
      pt := ptlist[ind];
      aa := pt[1];
      bb := pt[2];
      cc := pt[3];
      dd := pt[4];
      jnew := Evaluate(j,(aa*t+bb)/(cc*t+dd));
      changemat := changemat*Matrix([[aa,bb],[cc,dd]]);
      printf "After round %o size = %o.\n",it,minsize;
      j := jnew;
      //printf "j = %o.\n",j;
      if (it ge 5) then
        done := true;
      end if;
    else
      done := true;
    end if;
  end while;
  // Translation rounds
  printf "Doing two more random checks.\n";
  j, retM := transcheck(j);
  changemat := changemat*retM;
  // Hack for group number 7
  if size(Evaluate(j,t/31)) lt size(j) then
    j := Evaluate(j,t/31);
    changemat := changemat*Matrix([[1,0],[0,31]]);
  end if;
  printf "Final size is %o.\n",size(j);
  if size(j) ge prevsize then
    alldone := true;
    j := prevj;
    changemat := prevmat;
  else
    prevj := j;
    prevmat := changemat;
    prevsize := size(j);
  end if;
  end while;
  return j, changemat;
end function;

FF<tt> := FunctionField(Rationals());
RR<t> :=PolynomialRing(Rationals());
RRR<x,y> := PolynomialRing(Rationals(),2);
ZZ<b> := PolynomialRing(Integers());

// Compute the Weierstrass p-function to a precision of q^(prec/N).
// This is up to and includeing q^(prec/N).
// This routine returns the Fourier expansion
// of normalized version of p((c*t+d)/N;t)
// This is p multiplied by (9/Pi^2)
// We use the formula for the Fourier expansion from Chapter 4 of
// Diamond and Shurman. This is a holomorphic weight two modular
// form on Gamma(N). We look at ratios of these.

function weier(c,d,N,pr)
  c := c mod N;
  d := d mod N;
  const := R!(-3);
  if (c eq 0) then
    cosval := (1/2)*(zeta^d + zeta^(-d));
    const := const + R!(-18/(cosval - 1));
    pr := pr + 1;
  end if;
  ret := const;
  for m in [1..pr+N] do
    term := K!0;
    list1 := [ r*zeta^(d*r) : r in Divisors(m) | (Floor(m/r)-c) mod N eq 0 ];
    if #list1 gt 0 then
      term := term - 36*(&+list1);
    end if;
    list2 := [ r*zeta^(-d*r) : r in Divisors(m) | (Floor(m/r)-N+c) mod N eq 0];
    if #list2 gt 0 then
      term := term - 36*(&+list2);
    end if;
    if (m mod N eq 0) then
      term := term + 72*SumOfDivisors(Floor(m/N));
    end if;
    ret := ret + term*q^(m/N);
  end for;
  ret := ret + BigO(q^((pr+N+1)/N));
  return ret;
end function;

// This function takes two q-expansions and a degree bound
// The first one is a modular function, the second
// is a modular function of degree 1, and the
// third is a degree bound. This function tries to represent
// the first input as a rational function in the second.
// The two modular functions MUST belong to the same
// power series ring.
// The range of coefficients needed is from
// val2+deg*val1 to val2+deg*val1+2*deg+1

function ratfuncrep(modfunc,haup,deg)
  printf "Call to ratfuncrep with deg = %o.\n",deg;
  if IsWeaklyZero(modfunc) then
    return Parent(tt)!0;
  end if;
  field := BaseRing(Parent(modfunc));
  q := Parent(modfunc).1;
  M := ZeroMatrix(field,2*deg+2,2*deg+2);
  den := ExponentDenominator(haup);
  val1 := Min(0,Valuation(haup)*den);
  val2 := Min(0,Valuation(modfunc)*den);
  //printf "den = %o.\n",den;
  timet := Cputime();
  printf "Building last %o rows of the matrix.\n",deg+1;
  printf "Using coefficients %o to %o.\n",(val2+deg*val1)/den,(val2+deg*val1+2*deg+1)/den;
  for m in [0..deg] do
    haupprec := (m*val1 + 2*deg+2)/den;
    func2 := -(haup + BigO(q^(haupprec)))^(deg-m)*modfunc;
    //printf "For m = %o, the precision on func2 is from %o to %o.\n",m,Valuati\
on(func2),AbsolutePrecision(func2);
    //printf "For m = %o, precision needed is from %o to %o.\n",m,(val2+(deg-m)*val1)/den,(val2+(deg-m)*val1+2*deg+1)/den;
    //printf "Coefficient range %o to %o.\n",(val2+deg*val1)/den,(val2+deg*val1+2*deg+1)/den;
    for n in [1..2*deg+2] do
      M[m+deg+2][n] := Coefficient(func2,(val2+deg*val1+n-1)/den);
    end for;
  end for;
  printf "Time taken was %o.\n",Cputime(timet);
  timet := Cputime();
  printf "Building the first %o rows of the matrix.\n",deg+1;
  for m in [0..deg] do
    haupprec := (val2+m*val1+2*deg+2)/den;
    func2 := (haup+BigO(q^(haupprec)))^(deg-m);
    //printf "For m = %o, precision on func2 ranges from %o to %o.\n",m,Valuati\
on(func2),AbsolutePrecision(func2);
    //printf "Precision needed is %o to %o.\n",(val2+(deg-m)*val1)/den,(val2+(d\
eg-m)*val1+2*deg+1)/den;
    for n in [1..2*deg+2] do
      M[m+1][n] := Coefficient(func2,(val2+deg*val1+n-1)/den);
    end for;
  end for;
  printf "Time taken was %o.\n",Cputime(timet);
  printf "Solving the linear system of size %ox%o.\n",2*deg+2,2*deg+2;
  timet := Cputime();
  V := Nullspace(M);
  printf "Time taken was %o.\n",Cputime(timet);
  printf "Null space has dimension %o.\n",Dimension(V);
  v := V.1;
  QQ := Rationals();
  // We really hope that all the entries of V can be coerced into QQ
  numcoeffs := [ QQ!v[i] : i in [1..deg+1]];
  denomcoeffs := [ QQ!v[i] : i in [deg+2..2*deg+2]];
  mult := LCM([Denominator(v[i]) : i in [1..2*deg+2]]);
  numcoeffs := [ numcoeffs[i]*mult : i in [1..deg+1]];
  denomcoeffs := [ denomcoeffs[i]*mult : i in [1..deg+1]];
  ret := &+[ numcoeffs[i]*tt^(deg+1-i) : i in [1..deg+1]]/&+[ denomcoeffs[i]*tt\
^(deg+1-i) : i in [1..deg+1]];
  // New ret check
  numer := &+[ numcoeffs[i]*tt^(deg+1-i) : i in [1..deg+1]];
  denom := &+[ denomcoeffs[i]*tt^(deg+1-i) : i in [1..deg+1]];
  for d in [2..Dimension(V)] do
    vd := V.d;
    numerd := &+[ (QQ!vd[i])*tt^(deg+1-i) : i in [1..deg+1]];
    denomd := &+[ (QQ!vd[i])*tt^(2*deg+2-i) : i in [deg+2..2*deg+2]];
    if (numerd*denom - numer*denomd) ne 0 then
      printf "ERROR in ratfuncrep! is not unique!\n";
      bad := 0;
      bad2 := 1/bad;
    end if;    
  end for;
/*
  retcheck := &+[ numcoeffs[i]*(haup+BigO(q^AbsolutePrecision(modfunc)))^(deg+1-i) : i in [1..deg+1]]/&\
+[
  denomcoeffs[i]*(haup+BigO(q^AbsolutePrecision(modfunc)))^(deg+1-i) : i in [1..deg+1]];
  if IsWeaklyZero(retcheck - modfunc) then
    printf "Success! The linear system worked.\n";
    printf "The result was %o.\n",ret;
  else
    printf "Error! Solving the linear system failed!\n";
    bad := 0;
    bad2 := 1/bad;
  end if;
*/
  return ret;
end function;

// This function takes as input three modular functions
// modfunc, x, y. It also takes a degree bound deg.
// It then represents modfunc as a modular function
// of degree at most deg in the function field of the elliptic
// curve y^2 + a1 xy + a3 y = x^3 + a2 x^2 + a4 x + a6.

function genusoneratfuncrep(modfunc,x,y,deg)
  printf "Call to genusoneratfuncrep with deg = %o.\n",deg;
  field := BaseRing(Parent(modfunc));
  q := Parent(modfunc).1;
  den := LCM(ExponentDenominator(x),ExponentDenominator(y));
  valm := Valuation(modfunc)*den;
  valx := Valuation(x)*den;
  valy := Valuation(y)*den;
  // What range of exponents should we look at?
  val1 := valm + (3*deg-3)+Min(0,valx);
  val2 := (3*deg-3)*Min(0,valx);
  val3 := (3*deg-3)*Min(0,valx)+valy;
  startval := Min([val1,val2,val3]);
  rowlen := 13*deg;
  endval := startval + rowlen;
  //printf "The range of exponents is from %o to %o.\n",startval,endval;
  //printf "The actual precision needed for modfunc is %o.\n",(endval+Max(0,(3*deg-3)*(-valx)))/den;
  // An elliptic function of degree k has the form
  // (A(x) + B(x) y)/C(x), where deg A <= 3k-3, deg B <= 3k-5, and
  // deg C <= 3k-3.
  M := ZeroMatrix(field,9*deg-8,rowlen);
  printf "Building matrix - part 1.\n";
  for m in [0..3*deg-3] do
    relprec := 1+endval - (3*deg-3-m)*valx - valm;
    absprec := relprec - valx;
    func := -(x+BigO(q^(absprec/den)))^(3*deg-3-m)*modfunc;
    for n in [1..rowlen] do
      M[m+1][n] := Coefficient(func,(startval+n-1)/den);
    end for;
  end for;
  printf "Building matrix - part 2.\n";
  for m in [0..3*deg-3] do
    relprec := 1+endval - (3*deg-3-m)*valx;
    absprec := relprec - valx;
    func := (x+BigO(q^(absprec/den)))^(3*deg-3-m);
    for n in [1..rowlen] do
      M[3*deg-1+m][n] := Coefficient(func,(startval+n-1)/den);
    end for;
  end for;
  printf "Building matrix - part 3.\n";
  for m in [0..3*deg-6] do
    relprec := 1+endval - (3*deg-6-m)*valx - valy;
    absprec := relprec - valx;
    func := (x+BigO(q^(absprec/den)))^(3*deg-6-m)*y;
    for n in [1..rowlen] do
      M[6*deg-3+m][n] := Coefficient(func,(startval+n-1)/den);
    end for;
  end for;
  printf "Computing nullspace of %ox%o matrix.\n",NumberOfRows(M),NumberOfColumns(M);
  V := Nullspace(M);
  printf "Done!\n";
  printf "Null space has dimension %o.\n",Dimension(V);
  v := V.1;
  QQ := Rationals();
  // We really hope that all the entries of V can be coerced into QQ
/*
  for n in [1..9*deg-8] do
    printf "Solution %o is %o.\n",n,v[n];
  end for;
  bad := 0;
  bad2 := 1/bad;
*/
  FFF<ttt> := PolynomialRing(QQ);
  ccoeffs := [ QQ!v[i] : i in [1..3*deg-2]];
  acoeffs := [ QQ!v[i] : i in [3*deg-1..6*deg-4]];
  bcoeffs := [ QQ!v[i] : i in [6*deg-3..9*deg-8]];
  reta := &+[ acoeffs[i]*ttt^(3*deg-2-i) : i in [1..3*deg-2]];
  retb := &+[ bcoeffs[i]*ttt^(3*deg-5-i) : i in [1..3*deg-5]];
  retc := &+[ ccoeffs[i]*ttt^(3*deg-2-i) : i in [1..3*deg-2]];
  for d in [2..Dimension(V)] do
    vd := V.d;
    ccoeffsd := [ QQ!vd[i] : i in [1..3*deg-2]];
    acoeffsd := [ QQ!vd[i] : i in [3*deg-1..6*deg-4]];
    bcoeffsd := [ QQ!vd[i] : i in [6*deg-3..9*deg-8]];
    retad := &+[ acoeffsd[i]*ttt^(3*deg-2-i) : i in [1..3*deg-2]];
    retbd := &+[ bcoeffsd[i]*ttt^(3*deg-5-i) : i in [1..3*deg-5]];
    retcd := &+[ ccoeffsd[i]*ttt^(3*deg-2-i) : i in [1..3*deg-2]]; 
    chk1 := reta*retcd - retad*retc;
    chk2 := retb*retcd - retbd*retc;
    if (chk1 ne 0) or (chk2 ne 0) then
      printf "ERROR! Solution is not unique in genusoneratfuncrep.\n";
      bad := 0;
      bad2 := 1/bad;
    end if;
  end for;
  gcd := GCD([reta,retb,retc]);
  //printf "We found gcd = %o.\n",gcd;
  return reta div gcd, retb div gcd, retc div gcd;
end function;

// Function simplifyelliptic - This function takes reta, retb,
// and retc from the genusoneratfuncrep function, and also the
// ainvariants of the elliptic curve, and it outputs two polynomials
// p1 and p2 with the property that the input divided by p1(x) + p2(x)y
// is a square.

function simplifyelliptic(reta,retb,retc,ainvars)
  E := EllipticCurve(ainvars);
  F<x> := FunctionField(Rationals());
  S<zz> := PolynomialRing(F);
  f := zz^2 + ainvars[1]*zz*x + ainvars[3]*zz - (x^3 + ainvars[2]*x^2 + ainvars[4]*x + ainvars[5]);
  K<y> := FunctionField(f);
  elt := (K!Evaluate(reta,x)+K!Evaluate(retb,x)*y)/(K!Evaluate(retc,x));
  D := Divisor(elt);
  S := Support(D);
  infty := Divisor(InfinitePlaces(K)[1]);

  simpd := DivisorGroup(K)!0;
  D2 := DivisorGroup(K)!0;
  for i in [1..#S] do
    val := Valuation(D,S[i]);
    if (val mod 2 eq 1) then
      simpd := simpd + S[i];
    end if;
    D2 := D2 + Floor(val/2)*S[i];
  end for;
  deg := Degree(simpd);
  simpd := simpd - deg*infty;
  D2 := D2 - Degree(D2)*infty;

  // This divisor is probably not principal - Make it so it is!
  OK := EquationOrderFinite(K);

  // Guess a random point on E to add to the divisor to make it principal
  T, phi0 := TorsionSubgroup(E);
  torspoints := { phi0(a) : a in T };
  // We assume that E has rank 1.
  GG, phi := MordellWeilGroup(E);
  num := #Generators(GG);
  P := phi(GG.num);

  found := false;
  it := 0;
  while found eq false do
    chklist := [ it*P + s : s in torspoints ] cat [ -it*P + s : s in torspoints ];
    it2 := 0;
    done := false;
    while done eq false do
      it2 := it2 + 1;
      if it2 gt #chklist then
        done := true;
      else
        pt := chklist[it2];
        if (pt eq E!0) then
          dv := infty;
        else
          xcord := pt[1];
          ycord := pt[2];
          I := ideal<OK | x-xcord,y-ycord>;
          dv := Divisor(Place(I));
        end if;
        dv := dv - infty;
        if IsPrincipal(D2-dv) then
          done := true;
          found := true;
          //printf "We found it! We have pt = %o.\n",pt;
        end if;
      end if;
    end while;
    it := it + 1;
  end while;

  divg := D2-dv;
  a, gfunc := IsPrincipal(divg);
  g := Evaluate(RationalFunction(gfunc),y);
  simp1 := elt/g^2;
  T<z> := PolynomialRing(Rationals());
  elts := Eltseq(simp1);
  elts := [ T!elts[1], T!elts[2] ];
  if Degree(elts[1]) ge 0 then
    denomlcm1 := LCM([ Denominator(Coefficient(elts[1],i)) : i in [0..Degree(elts[1])]]);
    numgcd1 := GCD([ Numerator(Coefficient(elts[1],i)) : i in [0..Degree(elts[1])]]);
  else
    denomlcm1 := 1;
    numgcd1 := 0;
  end if;
  if Degree(elts[2]) ge 0 then
    denomlcm2 := LCM([ Denominator(Coefficient(elts[2],i)) : i in [0..Degree(elts[2])]]);
    numgcd2 := GCD([ Numerator(Coefficient(elts[2],i)) : i in [0..Degree(elts[2])]]);
  else
    denomlcm2 := 1;
    numgcd2 := 0;
  end if;
  denomlcm := LCM([denomlcm1,denomlcm2]);
  numgcd := GCD([numgcd1,numgcd2]);
  aa, bb := SquarefreeFactorization(denomlcm);
  cc, dd := SquarefreeFactorization(numgcd);  
  fac := (aa^2*bb^2)/dd^2;
  return elts[1]*fac,elts[2]*fac;
end function;

// This function act takes an element of the subspace of Eisenstein series
// represented by a vector giving the coefficients in terms of the
// basislist with entries in K, a level N, a matrix g and it computes the 
// image of this element under the matrix g.

// The action is on the right, as in Shimura.

function act(elt,g,N,basislist)
  newvec := [ K!0 : i in [1..#basislist]];
  det := lift(Determinant(g),N);
  for i in [1..#basislist-1] do
    a := basislist[i][1];
    b := basislist[i][2];
    imag := <lift(g[1][1]*a+g[2][1]*b,N),lift(g[1][2]*a+g[2][2]*b,N)>;
    if (imag[1] gt N/2) then
      imag[1] := N-imag[1];
      imag[2] := N-imag[2];
    end if;
    if ((imag[1] eq 0) or (imag[1] eq N/2)) and (imag[2] gt N/2) then
      imag[1] := N-imag[1];
      imag[2] := N-imag[2];
    end if;
    imag[1] := imag[1] mod N;
    imag[2] := imag[2] mod N;
    ind := Index(basislist,imag);
    if (ind eq 0) then
      printf "Error! The vector %o is not in our basis list!.\n",imag;
      bad := 0;
      bad2 := 1/bad; 
    end if;
    // Act on elt by zeta^i -> zeta^(i*det)
    coeffs := Eltseq(elt[i]);
    newvec[ind] := &+[ coeffs[i]*zeta^((i-1)*det) : i in [1..#coeffs]];
  end for;
  c := newvec[#basislist];
  for i in [1..#basislist] do
    newvec[i] := newvec[i] - c;
  end for;  
  newvec2 := [ newvec[i] : i in [1..#basislist-1]];
  return VectorSpace(K,#basislist-1)!newvec2;
end function;

// The function below uses the models and maps we have computed to
// map points down to the j-line. It takes two parameters -
// a subgroup number, and a list of rationals.

function maptoj(subnum,ptlist);
  curgp := subnum;
  curptlist := ptlist;
  while (curgp gt 1) do
    covergp := newsublist[curgp][7];
    mp := eval Read(Sprintf("%omap%o.txt",curgp,covergp));  
    curpt := Domain(mp)!curptlist;
    newpt := mp(curpt);
    curptlist := Coordinates(newpt);
    curgp := covergp;
  end while;
  return newpt;
end function;

// Step 2 - Find a nice subgroup to work with

if newsublist[cover][2] ne N then
  covergp := liftsub(newsublist[cover][3],Valuation(newsublist[cover][2],2),Valuation(N,2));
else
  covergp := newsublist[cover][3];
end if;
bigG := GL(2,Integers(N));
KK := sub<bigG | [3,0,0,1], [5,0,0,1], [1,1,0,1], [-1,0,0,-1]>;
phi, B, C := CosetAction(bigG,KK);
permcover := phi(covergp);

indbound := 8;
if subnum in [162,164,169,172,173,174,175,176,421,422,616,617,640,643,645,646,647,652,653,655,658,704,705,707,708,709,710,711,715,716,717] then
  indbound := 16;
end if;
if subnum in [177,178,441] then
  indbound := 24;
end if;
printf "Computing subgroups of the covering group of index <= %o.\n",indbound;
// These are only the subgroups up to conjugacy!!!
L0 := LowIndexSubgroups(permcover,indbound);
printf "Done! There were %o conjugacy classes.\n",#L0;

oldcuspnums := [];
printf "Computing number of cusps for each of these groups...";
for n in [1..#L0] do
  g, inde, cusps, e2, e3 := genus2(L0[n] @@ phi);
  Append(~oldcuspnums,cusps);  
end for;
printf "Done!\n";

printf "Building the list of all subgroups of index <= %o.\n",indbound;
L := [];
cuspnums := [];
for n in [1..#L0] do
  cursub := L0[n];
  Norm := Normalizer(permcover,cursub);
  cosets := RightTransversal(permcover,Norm);
  for m in [1..#cosets] do
    Append(~L,Conjugate(cursub,cosets[m]));
    Append(~cuspnums,oldcuspnums[n]);
  end for;
end for;
printf "Done. There are %o subgroups of index <= %o.\n",#L,indbound;

sub := newsublist[subnum][3];
permgp := phi(newsublist[subnum][3]);
printf "Finding which subgroups are contained in permgp...";
sublist := [ n : n in [1..#L] | L[n] subset permgp];
printf "Done!\n";

done := false;
found := false;
it := 1;
goodgroup := 0;
while (done eq false) do
  if (it in sublist) then
    //printf "Testing group %o.\n",it;
    cusps := cuspnums[it];
    llst := [ cuspnums[j] : j in [1..it-1] | L[it] subset L[j]];
    if #llst gt 0 then
      cuspsuper := Max(llst);
    else
      cuspsuper := 0;
    end if;
    if (cusps gt cuspsuper) then
      inde := Index(permgp,L[it]);
      sl2ind := Index(sub meet SL(2,Integers(N)),(L[it] @@ phi) meet SL(2,Integers(N)));
      if (inde eq sl2ind) then
        done := true;
        found := true;
        goodgroup := it;
        printf "Group %o works!\n",it;
        printf "It has %o cusps.\n",cusps;
      end if;
    end if;
    if (done eq false) then
      it := it + 1;
    end if;
  else
    it := it + 1;
  end if;
  if (it gt #L) then
    done := true;
    printf "We didn't find a good subgroup!\n";
    bad := 0;
    bad2 := 1/bad;
  end if;
end while;

// Here's the group we'll compute Eisenstein series for.
subsub := L[goodgroup] @@ phi;
permsubsub := L[goodgroup];

// Step 3 - Compute Eisenstein series

// First do the forms that transform nicely under SL_2.
// Note that there are 3*N^2/8 cusps on Gamma(N) if N is a power of 2, and so the space of Eisenstein
// series (which should be spanned by the xcoords) has dimension 3*N^2/8 - 1. The relation between these
// is that the sum of all of them is zero (since the sum of all of them is a holomorphic modular form of
// weight 2 on Gamma(N)).

basislist := [];
for a in [0..Floor(N/2)] do
  for b in [0..N-1] do
    if GCD([a,b,N]) eq 1 then
      addtolist := true;
      if (a eq 0) and (b gt N/2) then
        addtolist := false;
      end if;
      if (a eq N/2) and (b gt N/2) then
        addtolist := false;
      end if;
      if addtolist eq true then
        //printf "Adding %o to the basis list.\n",<a,b>;
        Append(~basislist,<a,b>);
      end if;
    end if;
  end for;
end for;
lastbasis := basislist[#basislist];
dimen := #basislist - 1;
printf "The space of Eisenstein series has dimension %o.\n",dimen;

// Represent an element of M_2(Gamma) that is a linear combination of the p_tau as a list
// of coefficients in Q(zeta) corresponding to the basis vectors.

printf "Computing action of generators of subsub intersect SL_2 on Eisenstein series for Gamma(%o).\n",N;

matlist := [];
sl2sub := subsub meet SL(2,Integers(N));
for g in Generators(sl2sub) do
  // Compute a matrix representing g on the space.
  M := ZeroMatrix(Rationals(),dimen,dimen);
  for i in [1..dimen] do
    // We could call act here, but it is wasteful since we only need the
    // image of a single basis element
    rowmat := ZeroMatrix(K,1,#basislist);
    a := basislist[i][1];
    b := basislist[i][2];
    imag := <lift(g[1][1]*a+g[2][1]*b,N),lift(g[1][2]*a+g[2][2]*b,N)>;
    if (imag[1] gt N/2) then
      imag[1] := N-imag[1];
      imag[2] := N-imag[2];
    end if;
    if ((imag[1] eq 0) or (imag[1] eq N/2)) and (imag[2] gt N/2) then
      imag[1] := N-imag[1];
      imag[2] := N-imag[2];
    end if;
    imag[1] := imag[1] mod N;
    imag[2] := imag[2] mod N;
    indind := Index(basislist,imag);
    if (indind eq 0) then
      printf "Error! The vector %o is not in our basis list!.\n",imag;
      bad := 0;
      bad2 := 1/bad; 
    end if;
    rowmat[1][indind] := 1;
    if rowmat[1][#basislist] ne 0 then
      c := rowmat[1][#basislist];
      for j in [1..#basislist] do
        rowmat[1][j] := rowmat[1][j] - c;
      end for;
    end if;
    for j in [1..dimen] do
      M[j][i] := rowmat[1][j];
    end for;
  end for;
  Append(~matlist,M);
end for;

I := IdentityMatrix(Rationals(),dimen);
if #matlist eq 0 then
  Append(~matlist,I);
end if;
printf "Finding vectors invariant under subsub intersect SL_2.\n",subnum;
Ker := Kernel(Transpose(matlist[1]-I));
for m in [2..#matlist] do
  Ker := Ker meet Kernel(Transpose(matlist[m]-I));
end for;
kerbasis := Basis(Ker);
printf "The dimension of the space of Eisenstein series for subsub intersect SL_2 is %o.\n",#kerbasis;
if (#kerbasis ne (cusps - 1)) then
  printf "Error! We didn't get enough Eisenstein series!.\n";
end if;

// Now we look at the Q-span of the forms with coefficients in Q(zeta_N)
// that are linear combinations of the elements in K.

phiN := EulerPhi(N);
printf "Computing the action of the generators of subsub on the Q-vector space of dimension %o.\n",phiN*#kerbasis;
matlist := [];
genlist := SetToSequence(Generators(subsub));
for g in genlist do
  M := ZeroMatrix(Rationals(),phiN*#kerbasis,phiN*#kerbasis);
  for k in [1..#kerbasis] do
    // Compute image of zeta^i*(kth basis vector) under the action of g
    basisvec := act(kerbasis[k],g,N,basislist);
    //printf "Basis vector is %o.\n",[ basisvec[j] : j in [1..dimen]];
    //printf "Attempting to coerce into K.\n";
    kelt := Ker![ basisvec[j] : j in [1..dimen]];
    coords := Coordinates(Ker,kelt);

    det := Integers()!Determinant(g);
    for i in [0..phiN-1] do
      // Write the image of zeta^i*(kth basis vector) to column
      // (i+1) + phiN*k
      fac := Eltseq(zeta^(i*det));
      for j in [0..phiN-1] do
        for m in [1..#kerbasis] do
          entry := coords[m]*fac[j+1];
          M[(m-1)*phiN+(j+1)][i+1 + phiN*(k-1)] := entry;
        end for;
      end for;
    end for;
  end for;
  Append(~matlist,M);
end for;

I := IdentityMatrix(Rationals(),phiN*#kerbasis);
if #matlist eq 0 then
  Append(~matlist,I);
end if;
printf "Finding vectors invariant under subsub.\n",subnum;
Ker2 := Kernel(Transpose(matlist[1]-I));
for m in [2..#matlist] do
  Ker2 := Ker2 meet Kernel(Transpose(matlist[m]-I));
end for;
ker2basis := Basis(Ker2);
printf "The dimension is %o.\n",#ker2basis;

// Translate back to the original form

finalbasismat := ZeroMatrix(K,#ker2basis,dimen);
Vsp := VectorSpace(K,dimen);
for n in [1..#ker2basis] do
  vec := Vsp!0;
  for m in [1..phiN*#kerbasis] do
    formnum := Floor((m-1)/phiN) + 1;
    rootnum := m - phiN*(formnum-1) - 1;
    vec := vec + ker2basis[n][m]*zeta^rootnum*(Vsp!kerbasis[formnum]);
  end for;
  finalbasismat[n] := vec;
end for;

// Enumerate the representatives of the right cosets of subsub in covergp 
// sorted by coset of sub in covergp

// S1 and S2 are sets of representatives for the right cosets

S1 := SetToSequence(RightTransversal(covergp,sub));
S2 := SetToSequence(RightTransversal(sub,subsub));

bigind := Index(covergp,subsub);
printf "The index of subsub in covergp is %o.\n",bigind;

seq := [ i : i in [1..#ker2basis]];
formtry := &+[ seq[m]*finalbasismat[m] : m in [1..#ker2basis]];

printf "We have selected a 'random' element of the Eisenstein subspace.\n";

imagelist := [];
images := {};
for i in [1..#S1] do
  imlist := [];
  for j in [1..#S2] do
    im := act(formtry,S2[j]*S1[i],N,basislist);
    Include(~images,im);
    Append(~imlist,im);
  end for;
  Append(~imagelist,imlist);
end for;
printf "It has %o images under the action of covergp.\n",#images;
if #images ne bigind then
  printf "This is a serious problem.\n";
  printf "Error! It doesn't have enough images!.\n";
  bad := 0;
  bad2 := 1/bad;
end if;

// Step 4 - Compute Fourier expansions

printf "Computing Weierstrass p-function Fourier expansions.\n";
xcoords := ZeroMatrix(R,N,N);
for i in [1..#basislist-1] do
  a := basislist[i][1];
  b := basislist[i][2];
  //printf "Computing expansion %o of %o.\n",i,#basislist-1;
  xcoords[a+1][b+1] := weier(a,b,N,prec);
end for;

fourierlist := [];
fouriernum := 0;
for i in [1..#S1] do
  ilist := [];
  for j in [1..#S2] do
    fouriernum := fouriernum + 1;
    printf "Computing Fourier expansion %o.\n",fouriernum;
    fourier := &+[ imagelist[i][j][k]*xcoords[basislist[k][1]+1][basislist[k][2]+1] : k in [1..#basislist-1]];
    Append(~ilist,fourier);
  end for;
  Append(~fourierlist,ilist);
end for;

printf "Symmetrizing.\n";
wt := 0;
formsused := [];
if (subsub eq sub) then
  wt := 2;
  formsused := [ fourierlist[i][1] : i in [1..#fourierlist]];
else
  wt := 4;
  termnum := 0;
  for i in [1..#S1] do
    formtouse := R!0;
    maxnum := #S1*#S2*(#S2-1)/2;
    for j in [1..#S2] do
      for j2 in [j+1..#S2] do
        termnum := termnum+1;
        printf "Doing term %o of %o.\n",termnum,maxnum;
        formtouse := formtouse + fourierlist[i][j]*fourierlist[i][j2];
      end for;
    end for;
    Append(~formsused,formtouse);
  end for;
end if;
printf "Done!\n";
chk := #{ formsused[i] : i in [1..#formsused]};
if (chk lt #S1) then
  printf "Error - the number of roots is %o. This is less than %o.\n",chk,#S1;
  printf "Trying again.\n";
  wt := 6;
  formsused := [];
  termnum := 0;
  for i in [1..#S1] do
    formtouse := R!0;
    maxnum := #S1*(#S2)*(#S2-1)*(#S2-2)/6;
    for j in [1..#S2] do
      for j2 in [j+1..#S2] do
        for j3 in [j2+1..#S2] do
          termnum := termnum + 1;
          printf "Doing term %o of %o.\n",termnum,maxnum;
          formtouse := formtouse + fourierlist[i][j]*fourierlist[i][j2]*fourierlist[i][j3];
        end for;
      end for;
    end for;
    Append(~formsused,formtouse);
  end for;
  rechk := #{ formsused[i] : i in [1..#formsused]};
  printf "This time we have %o roots.\n",rechk;
  if (rechk lt #S1) then
    printf "We're still screwed!\n";
    bad := 0;
    bad2 := 1/bad;
  end if;
end if;

effwt := 0;
if (wt eq 6) then
  effwt := 6;
  denomfunc := Eisenstein(6,q : Precision := Ceiling(2*prec/N));
end if;
if (wt eq 4) then
  effwt := 4;
  denomfunc := Eisenstein(4,q : Precision := Ceiling(2*prec/N));
end if;
if (wt eq 2) then
  if (newsublist[cover][1] mod 3 eq 0) then
    effwt := 2;   
    denomfunc := 2*Eisenstein(2,q^2 : Precision := Ceiling(2*prec/N)) - Eisenstein(2,q : Precision := Ceiling(2*prec/N));
  else
    effwt := 6;
    denomfunc := Eisenstein(6,q : Precision := Ceiling(2*prec/N))/Eisenstein(4,q : Precision := Ceiling(2*prec/N));
  end if;    
end if;

// Step 5 - Compute the minimal polynomial over Q(X(cover)) of
// the modular function formsused[1]/denomfunc.

degbound := Ceiling((effwt*newsublist[subnum][1])/(12*ind));
S<xxx> := PolynomialRing(R);
poly := S!1;
for i in [1..#S1] do
  poly := poly*(denomfunc*xxx - formsused[i]);
end for;
poly := poly/(denomfunc^Degree(poly));
modf := formsused[1]/denomfunc;

// Step 6 - Now that we have a model, work with it!
// What we do now depends on the genus of the cover

// Here we initialize two flags that will control what we output. The first
// will say if the modular curve has infinitely many rational points or not.

inftyratpoints := false;

// The second flag will be indicate if there are no rational points on the
// curve or not. If we can prove that there are no points, then we don't
// need to output the model and the covering map, otherwise we do. 

noratpoints := false;

gen := newsublist[subnum][5];

if gencover eq 0 then
  coefflist := [];
  for m in [0..Degree(poly)-1] do
    printf "Call to ratfunc with m = %o.\n",m;
    ret := ratfuncrep(Coefficient(poly,m),haupt,(Degree(poly)-m)*degbound);
    //printf "Return of ratfunc is %o.\n",ret;
    Append(~coefflist,ret);
    printf "Result has degree %o.\n",Max(Degree(Numerator(ret)),Degree(Denominator(ret)));
  end for;
  Append(~coefflist,1);
  denom := FF!LCM([Denominator(coefflist[i]) : i in [1..#coefflist]]);
  scaledcoefflist := [ RR!(denom*coefflist[i]) : i in [1..#coefflist]];
  curveeq := &+[ Evaluate(scaledcoefflist[i+1],x)*y^i : i in [0..#coefflist-1]];
  cofs := Coefficients(curveeq,y);
  if not IsIrreducible(curveeq) then
    printf "Error! We didn't get a good primitive element!\n";
    printf "The factorization is %o.\n",Factorization(curveeq);
    bad := 0;
    bad2 := 1/bad;
  end if;
  // Case 1 - Degree is 2
  if (ind eq 2) then
    // First simplify the square root
    cofs := Coefficients(curveeq,y);
    disc := RR!Evaluate(cofs[2]^2 - 4*cofs[1]*cofs[3],[t,0]);
    denom := LCM([Denominator(Coefficient(disc,m)) : m in [0..Degree(disc)]]);
    FFF := Factorization(denom);
    if #FFF gt 0 then
      mult := &*[ FFF[n][1]^(2*Ceiling(FFF[n][2]/2)) : n in [1..#FFF]];
      disc := disc*mult;
      disc := ZZ!disc;
    else
      disc := ZZ!disc;
    end if;
    FFF2 := Factorization(disc);
    if #FFF2 gt 0 then
      discdiv := &*[ FFF2[n][1]^(2*Floor(FFF2[n][2]/2)) : n in [1..#FFF2]];
      newdisc := disc div discdiv;
    else
      newdisc := disc;
    end if;
    printf "Fairly simple model is y^2 = %o.\n",newdisc;
    // Here, y is the new uniformizer, and newdisc is a polynomial in
    // the old one.
    // Case 1a - Degree is 2 and genus is zero
    if (gen eq 0) then
      // uni is the modular function equal to y.
      // We should be able to manually build a conic now
      a := Rationals()!Coefficient(newdisc,2);
      b := -Rationals()!1;
      c := Rationals()!Coefficient(newdisc,0);
      d := Rationals()!0;
      e := Rationals()!0;
      f := Rationals()!Coefficient(newdisc,1);
      CC := Conic([a,b,c,d,e,f]);
      isrational, P := HasRationalPoint(CC);
      printf "The modular curve is a conic given by %o.\n",DefiningEquation(CC);
      if isrational then
        inftyratpoints := true;
        printf "This genus zero curve has a rational point on it, namely %o.\n",P;
        A := Parametrization(CC,P);
        xparam := A([tt,1])[1];
        yparam := A([tt,1])[2];
        printf "We have x = %o.\n",xparam;
        printf "We have y = %o.\n",yparam;
        // Build a polynomial ring over FF, the rational function field.
        RRRR<xx> := PolynomialRing(FF);
        newcurveeq := -xx^2 + a*tt^2 + f*tt + c;
        SSSS<yy> := FunctionField(newcurveeq);
        TTTT<zz> := PolynomialRing(SSSS);
        num := Numerator(yparam);
        denom := Denominator(yparam);
        uniformpoly := Evaluate(denom,zz)*yy - Evaluate(num,zz);
        unifunc := Roots(uniformpoly)[1][1];

        // unifunc is a uniformizer for the function field
        newpoly := MinimalPolynomial(unifunc);
        denom := LCM([Denominator(Coefficient(newpoly,i)) : i in [0..Degree(newpoly)]]);
        newpoly2 := &+[ Evaluate(RR!(denom*Coefficient(newpoly,i)),y)*x^i : i in [0..Degree(newpoly)]];
        if not Degree(newpoly2,y) eq 1 then
          printf "Something is wrong! We don't have a uniformizer!\n";
          printf "The polynomial is %o.\n",newpoly2;
          bad := 0;
          bad2 := 1/bad;
        end if;
        newcofs := Coefficients(newpoly2,y);
        initratfunc := -Evaluate(newcofs[1],[tt,0])/Evaluate(newcofs[2],[tt,0]);
        printf "We have y = %o.\n",initratfunc;
        printf "Here tt is the Hauptmodul for the group %o,\n",subnum;
        printf "and y is the Hauptmodul for group %o.\n",cover;
        printf "Optimizing.\n";
        newratfunc, M := supersimplify(initratfunc);
        printf "The final model is y = %o.\n",newratfunc;
        M := ChangeRing(M,Rationals());
        M := M^(-1);
        if subnum in covergpset then
          // Tweak uniformizer
          printf "Computing modular function.\n";
          truncprec := Ceiling(AbsolutePrecision(modf)+20);
          // cofs := Coefficients(curveeq,y);
          modfunccofs := [ Evaluate(RR!Evaluate(cofs[i],[t,0]),(haupt+BigO(q^truncprec))): i in [1..#cofs]];
          olduni := 2*modfunccofs[3]*modf+modfunccofs[2];
          if #FFF gt 0 then
            aa, multsqrt := IsSquare(mult);
            olduni := olduni*multsqrt;
          end if;
          if #FFF2 gt 0 then
            aa, discdivsqrt := IsSquare(discdiv);
            uni := olduni/Evaluate(discdivsqrt,haupt+BigO(q^truncprec));
          else
            uni := olduni;
          end if;
          seqlist := ElementToSequence(unifunc);
          truncprec := Ceiling(AbsolutePrecision(uni)+20);
          printf "Computing uniformizer coefficients.\n";
          coefflist := [ Evaluate(seqlist[i],haupt+BigO(q^truncprec)) : i in [1..#seqlist]];
          printf "Computing uniformizer...\n";
          uniformize := &+[ coefflist[i]*uni^(i-1) : i in [1..#seqlist]];
          newhaupt := (M[1][1]*uniformize + M[1][2])/(M[2][1]*uniformize + M[2][2]);
        end if;
      else
        printf "This genus zero curve has no rational points!\n";
        noratpoints := true;  
      end if;
    end if;
    // Case 1b - Degree is 2 and genus is one
    if (gen eq 1) then
      locsov := true;
      if Degree(newdisc) eq 4 then
        HypC := HyperellipticCurve(newdisc);
        printf "Analyzing the curve y^2 = %o.\n",newdisc;
        if IsLocallySolvable(HypC,2) then
          printf "The curve has 2-adic points.\n";
          if (QuarticNumberOfRealRoots(newdisc) ge 1) or (LeadingCoefficient(newdisc) gt 0) then
            printf "The curve has real points too.\n";
            printf "Therefore it probably has rational points.\n";
            // All the modular curves we've computed have good reduction away
            // from 2. Note that the Shafarevich-Tate groups of all elliptic
            // curves with good reduction away from 2 over Q are trivial,
            // at least assuming strong BSD.
          else
            printf "The curve has no real points!\n";
            locsov := false;
          end if;
        else
          printf "The curve does not have any 2-adic points, and so it doesn't have any rational points either!\n";
          locsov := false;
        end if;
        if (locsov eq false) then
          noratpoints := true;
        else
          // If we get here, the curve should have a rational point.
          pointfound := false;
          searchbound := 1;
          pt := 0;
          while pointfound eq false do
            searchbound := searchbound * 10;
            S := Points(HypC : Bound := searchbound);
            if #S gt 0 then
              pointfound := true;
              pt := SetToIndexedSet(S)[1];
            end if;
          end while;
          printf "We found a point. It is %o.\n",pt;
          E, phi0 := EllipticCurve(HypC,pt);
          s := DefiningPolynomials(phi0);
          E2, phi := MinimalModel(E);
          printf "The modular curve is elliptic. It is the %o.\n",E2;
          s2 := DefiningPolynomials(Inverse(phi0*phi));
        end if;
      end if;
      if Degree(newdisc) eq 3 then
        newdisc := Evaluate(newdisc,t);
        // Right now, x is the old Hauptmodul, and y is the generator
        // of the new
        HypC := HyperellipticCurve(newdisc);
        E, phi0 := EllipticCurve(HypC);
        s := DefiningPolynomials(phi0);
        E2, phi := MinimalModel(E);
        printf "The modular curve is elliptic. It is the %o.\n",E2;
        s2 := DefiningPolynomials(Inverse(phi0*phi));
      end if;
      // Some writing data to files was removed here
      if (locsov eq true) then
        GG, phiphi := MordellWeilGroup(E2);
        if AnalyticRank(E2) eq 0 then
          printf "The rank of this elliptic curve is zero and the rational points are the following.\n";
          // Create the map from Xcur -> Xcover
          EE<XX,YY,ZZ> := EllipticCurve(aInvariants(E2));
          poly1 := Evaluate(s2[1],[XX,YY,ZZ]);
          poly2 := Evaluate(s2[3],[XX,YY,ZZ]);
          Xcov := eval Read("eq" cat IntegerToString(cover) cat ".txt");
          maptocover := map<EE -> Xcov | [poly1,poly2]>;
          for a in GG do
            curpt := phiphi(a);
            jimage := maptoj(cover,Coordinates(maptocover(curpt)));
            if jimage[2] eq 0 then
              printf "The point %o is a cusp.\n",curpt;
            else
              printf "The point %o maps to j = %o.\n",curpt,jimage[1];
            end if;
          end for;
        else
          inftyratpoints := true;
          printf "The rank of this elliptic curve is positive.\n";
          printf "The Mordell-Weil group is %o.\n",GG;
          printf "Generators are:\n";
          for a in Generators(GG) do
            printf "%o.\n",phiphi(a);
          end for;
          if subnum in covergpset then
            printf "Computing modular function.\n";
            truncprec := Ceiling(AbsolutePrecision(modf)+20);
            modfunccofs := [ Evaluate(RR!Evaluate(cofs[i],[t,0]),(haupt+BigO(q^truncprec))): i in [1..#cofs]];
            olduni := 2*modfunccofs[3]*modf+modfunccofs[2];
            if #FFF gt 0 then
              aa, multsqrt := IsSquare(mult);
              olduni := olduni*multsqrt;
            end if;
            if #FFF2 gt 0 then
              aa, discdivsqrt := IsSquare(discdiv);
              uni := olduni/Evaluate(discdivsqrt,haupt+BigO(q^truncprec));
            else
              uni := olduni;
            end if;
            xmodfunc := haupt;
            ymodfunc := uni;
            firstx := Evaluate(s[1],[xmodfunc,ymodfunc,1]);
            firsty := Evaluate(s[2],[xmodfunc,ymodfunc,1]);
            firstz := Evaluate(s[3],[xmodfunc,ymodfunc,1]);
            xmodfunc := firstx/firstz;
            ymodfunc := firsty/firstz;
            s3 := DefiningPolynomials(phi);
            nextx := Evaluate(s3[1],[xmodfunc,ymodfunc,1]);
            nexty := Evaluate(s3[2],[xmodfunc,ymodfunc,1]);
            nextz := Evaluate(s3[3],[xmodfunc,ymodfunc,1]);
            xmodfunc := nextx/nextz;
            ymodfunc := nexty/nextz;
          end if;
        end if;
      end if;
    end if;
    // Case 1c - Degree is 2 and genus is >= 2 (only genus 2 happens here)
    if (gen ge 2) then
      HypC := HyperellipticCurve(newdisc);
      newC, phi := ReducedMinimalWeierstrassModel(HypC);
      printf "Reduced model is %o.\n",newC;
      printf "Map is %o.\n",DefiningPolynomials(phi);
      J := Jacobian(HypC);
      highbound := RankBound(J);
      printf "The rank of the Jacobian is at most %o.\n",highbound;
      if (highbound eq 0) then
        PS := Chabauty0(J);
        PS := { phi(pt) : pt in PS };
        printf "The following is a list of all the points on the curve:\n";
      else
        PS := Points(newC : Bound := 100000);
        printf "Here are some points on the curve that we found:\n";
      end if;
      // Build the map to the cover
      phi2 := Inverse(phi);
      for pt in PS do
        pt2 := phi2(pt);
        printf "The point on the cover is %o.\n",pt2;
        if (pt2[3] eq 0) then
          p1pt := [1,0]; 
        else
          p1pt := [pt2[1]/pt2[3],1];
        end if;
        jinvar := maptoj(newsublist[subnum][7],p1pt);
        if (jinvar[2] eq 0) then
          printf "The point %o is a cusp.\n",pt;
        else
          printf "The point %o maps to j = %o.\n",pt,jinvar[1];
        end if;
        if IsSingular(newC!pt) then
          printf "The latter point is singular.\n";
        end if;
      end for;
    end if;
  end if;
  // Case 2 - Degree is > 2
  if (ind gt 2) then
    C := ProjectiveClosure(Curve(AffineSpace(Rationals(),2),curveeq));
    if (gen eq 0) then
      // This is called only for X_6, X_7 and X_55
      bound := 10;
      found := false;
      while found eq false do
        printf "Searching for a non-singular rational point with bound %o.\n",bound;
        plist := PointSearch(C,bound);
        it := 1;
        done := false;
        while done eq false do
          pt := plist[it];
          if not IsSingular(C!pt) then
            P := C!pt;
            done := true;
            found := true;
          else
            it := it + 1;
          end if;
          if (it gt #plist) then
            done := true;
            bound := bound*2;
          end if;
        end while;
        if bound gt 60000 then
          printf "No point found. Giving up.\n";
          bad := 0;
          bad2 := 1/bad;
        end if;
      end while;
      inftyratpoints := true;
      printf "Parametrizing the curve.\n";
      A := Parametrization(C,P);
      printf "Found it!\n";
      xparam := A([tt,1])[1];
      yparam := A([tt,1])[2];
      printf "We have x = %o.\n",xparam;
      printf "We have y = %o.\n",yparam;
      // Build a polynomial ring over FF, the rational function field.
      RRRR<xx> := PolynomialRing(FF);
      SSSS<yy> := FunctionField(Evaluate(curveeq,[tt,xx]));
      TTTT<zz> := PolynomialRing(SSSS);
      num := Numerator(yparam);
      denom := Denominator(yparam);
      uniformpoly := Evaluate(denom,zz)*yy - Evaluate(num,zz);
      unifunc := Roots(uniformpoly)[1][1];
      seqlist := ElementToSequence(unifunc);
      truncprec := Ceiling(AbsolutePrecision(modf)+20);
      printf "Computing uniformizer coefficients.\n";
      coefflist := [ Evaluate(seqlist[i],haupt+BigO(q^truncprec)) : i in [1..#seqlist]];
      printf "Computing uniformizer...\n";
      uniformize := &+[ coefflist[i]*modf^(i-1) : i in [1..#seqlist]];
      newpoly := MinimalPolynomial(unifunc);
      denom := LCM([Denominator(Coefficient(newpoly,i)) : i in [0..Degree(newpoly)]]);
      newpoly2 := &+[ Evaluate(RR!(denom*Coefficient(newpoly,i)),y)*x^i : i in [0..Degree(newpoly)]];
      if not Degree(newpoly2,y) eq 1 then
        printf "Something is wrong! We don't have a uniformizer!\n";
        printf "The polynomial is %o.\n",newpoly2;
        bad := 0;
        bad2 := 1/bad;
      end if;
      cofs := Coefficients(newpoly2,y);
      initratfunc := -Evaluate(cofs[1],[tt,0])/Evaluate(cofs[2],[tt,0]);
      printf "We have y = %o.\n",initratfunc;
      printf "Here tt is the Hauptmodul for the group %o,\n",subnum;
      printf "and y is the Hauptmodul for group %o.\n",cover;
      printf "Optimizing.\n";
      newratfunc, M := supersimplify(initratfunc);
      printf "The final model is y = %o.\n",newratfunc;
      //printf "M = %o.\n",M;
      M := ChangeRing(M,Rationals());
      M := M^(-1);
      newhaupt := (M[1][1]*uniformize + M[1][2])/(M[2][1]*uniformize + M[2][2]);
    end if;
    if (gen gt 0) then
      // This only happens for X_441 - we'll see if we can do the computation
      // anyway.
      printf "We must be doing X_441. Burcu Baran already worked this case out.\n";
      printf "The curve equation is %o.\n",C;
      printf "Determining hyperelliptic model.\n";
      A, B, phi0 := IsHyperelliptic(C);
      printf "The hyperelliptic model is %o.\n",B;
      printf "The map is %o.\n",phi0;
      // The answer should be yes.
      newC, phi := ReducedMinimalWeierstrassModel(B);
      printf "The minimal model is %o.\n",newC;
      printf "The map to the minimal model is %o.\n",phi;
      printf "The covering map written to a file will be wrong!\n";
    end if;  
  end if;
end if;

if gencover eq 1 then
  // In this case ind = 2
  // There are cases where gen = 1 here, but all of the base curves
  // have rank zero. That's why inftyratpoints = false.
  modfuncdisc := Coefficient(poly,1)^2 - 4*Coefficient(poly,0);
  reta, retb, retc := genusoneratfuncrep(modfuncdisc,xcoord,ycoord,2*degbound);
  newreta, newretb := simplifyelliptic(reta,retb,retc,alist);
  // OK, now we build the curve as a cover of the elliptic curve.
  deg := Max([2, 1 + Degree(newretb), Degree(newreta)]);
  P3 := ProjectiveSpace(Rationals(),3);
  CR<x,y,z,w> := CoordinateRing(P3);
  // This is where the curve should go
  // here we refer to alist, the a-invarints read from the model
  // of the cover.
  f1 := y^2*z + alist[1]*x*y*z + alist[3]*y*z^2
   - (x^3 + alist[2]*x^2*z + alist[4]*x*z^2 + alist[5]*z^3);
  f2 := w^2*z^(deg-2);
  if (newreta ne 0) then
    f2 := f2 - &+[ Coefficient(newreta,i)*x^i*z^(deg-i) : i in [0..Degree(newreta)]];
  end if;
  if (newretb ne 0) then
    f2 := f2 - &+[ Coefficient(newretb,i)*y*x^i*z^(deg-1-i) : i in [0..Degree(newretb)]];
  end if;
  printf "Our model is given by:\n";
  printf "%o = 0.\n",Evaluate(f1,[x,y,1,0]);
  printf "w^2 = %o + (%o)*y\n",Evaluate(newreta,x),Evaluate(newretb,x);
  I := ideal<CR | f1,f2>;
  I2 := PrimaryDecomposition(I);
  indextouse := Index([ IsPrime(I2[i]) : i in [1..#I2]],true);
  //printf "The primary decomposition is %o.\n",I2; 
  // Hopefully this is the right piece to look at.
  C := Curve(P3,I2[indextouse]); 
  printf "Computing genus of the model we got.\n"; 
  gencomp := Genus(C);
  printf "It is %o.\n",gencomp;
  if (gen ne gencomp) then
    printf "Error! Our model computation is wrong!.\n";
    bad := 0;
    bad2 := 1/bad;
  end if;
  // Do a search for rational points using the Mordell-Weil group of E
  mwbound := 100;
  EEE := EllipticCurve(alist);
  T, maptotors := TorsionSubgroup(EEE);
  torspoints := { maptotors(t) : t in T };
  MW, maptoMW := MordellWeilGroup(EEE);
  gener := maptoMW(MW.NumberOfGenerators(MW));
  // We are using that all of our elliptic curves have rank 1.
  bigptlist := [ m*gener + t : m in [-mwbound..mwbound], t in torspoints ];
  foundpoints := [];
  for p in bigptlist do
    if p[3] ne 0 then
      wsqr := Evaluate(newreta,p[1])+p[2]*Evaluate(newretb,p[1]);
      issqr, sqr := IsSquare(wsqr);
      if issqr then
        pt := [p[1],p[2],p[3],sqr];
        chk1, cpoint := IsCoercible(C,pt);
        if (chk1 eq true) then
          Append(~foundpoints,cpoint);
        end if;
        if (sqr ne 0) then
          pt := [p[1],p[2],p[3],-sqr];
          chk1, cpoint := IsCoercible(C,pt);
          if (chk1 eq true) then
            Append(~foundpoints,cpoint);
          end if;
        end if;     
      end if;
    end if;
  end for; 
  smallsearch := PointSearch(C,100);
  for pt in smallsearch do
    cpoint := C!pt;
    if Index(foundpoints,cpoint) eq 0 then
      Append(~foundpoints,cpoint);
    end if;
  end for; 
  printf "Here are the points we found on C:\n";
  for pt in foundpoints do
    if (pt[1] ne 0) or (pt[2] ne 0) or (pt[3] ne 0) then
      jinvar := maptoj(newsublist[subnum][7],[pt[1],pt[2],pt[3]]);
      if (jinvar[2] eq 0) then
        printf "The point %o is a cusp.\n",pt;
      else
        printf "The point %o maps to j = %o.\n",pt,jinvar[1];
      end if;
      if IsSingular(pt) then
        printf "The latter point is singular.\n";
      end if;
    else 
      printf "The point %o doesn't map to the elliptic curve.\n",pt;
      if IsSingular(pt) then
        printf "It is singular.\n";
      end if;
    end if;
  end for;
  if (gen eq 3) then
  // See if it is hyperelliptic or get the plane quartic model
  // but don't do anything with it.
    B := 0;
    A, B, phi2 := IsGeometricallyHyperelliptic(C);
    if (A eq true) then
      if HasRationalPoint(B) then
        printf "The curve is hyperelliptic. A hyperelliptic model is %o.\n",B;
        printf "The map is %o.\n",phi2;
      else
        printf "This genus 3 curve has a degree 2 map to a conic with no rational points!\n";
        noratpoints := true;
      end if;
    end if;
    if (A eq false) then
      printf "The curve is not geometrically hyperelliptic.\n";
      phi3 := CanonicalMap(C);
      C2 := CanonicalImage(C,phi3);
      quartmod := DefiningEquations(C2)[1];
      printf "It's plane quartic model is %o.\n",Evaluate(quartmod,[x,y,z]);
    end if;
  end if;  
end if;

// Step 7 - Do all the file writing once at the end.

printf "Noratpoints = %o.\n",noratpoints;
if (noratpoints eq false) then
  // Step 7a - Write the model
  filename := "eq" cat IntegerToString(subnum) cat ".txt"; 
  System("rm " cat filename);
  if (gen eq 0) then
    PrintFile(filename,Sprintf("X%o<a%o,b%o> := ProjectiveSpace(Rationals(),1);",subnum,subnum,subnum));
  end if;
  if (gen eq 1) then
    PrintFile(filename,Sprintf("X%o<x%o,y%o,z%o> := EllipticCurve(%o);",subnum,subnum,subnum,subnum,aInvariants(E2)));
  end if;
  if (gen eq 2) then
    f, h := HyperellipticPolynomials(newC);
    ABCDEF<l> := PolynomialRing(Rationals());
    f := ABCDEF!f;
    h := ABCDEF!h;
    PrintFile(filename,"ABCDEF<l> := PolynomialRing(Rationals());");
    PrintFile(filename,Sprintf("X%o<x%o,y%o,z%o> := HyperellipticCurve([%o,%o]);",subnum,subnum,subnum,subnum,f,h));   
  end if; 
  if (gen ge 3) then
    ABCDEF<X,Y,Z,W> := PolynomialRing(Rationals(),4);
    lst := DefiningPolynomials(C);
    lst := [ ABCDEF!lst[i] : i in [1..#lst]];
    PrintFile(filename,"P3<X,Y,Z,W> := ProjectiveSpace(Rationals(),3);");
    PrintFile(filename,Sprintf("X%o<x%o,y%o,z%o,w%o> := Curve(P3,",subnum,subnum,subnum,subnum,subnum));
    PrintFile(filename,Sprintf("%o);",lst));
  end if;
  PrintFile(filename,Sprintf("return X%o;",subnum));

  // Step 7b - Write the covering map

  filename := IntegerToString(subnum) cat "map" cat IntegerToString(cover)
  cat ".txt";
  System("rm " cat filename);
  Xsub := eval Read("eq" cat IntegerToString(subnum) cat ".txt");
  PrintFile(filename,Sprintf("X%o := eval Read(\"%o.txt\");",cover,"eq" cat IntegerToString(cover)));
  PrintFile(filename,Sprintf("X%o := eval Read(\"%o.txt\");",subnum,"eq" cat IntegerToString(subnum)));
  if (gen eq 0) then
    var1 := Xsub.1;
    var2 := Xsub.2;
    PrintFile(filename,Sprintf("%o := X%o.1;",var1,subnum));
    PrintFile(filename,Sprintf("%o := X%o.2;",var2,subnum));
    num := Numerator(newratfunc);
    denom := Denominator(newratfunc);
    degdeg := Max(Degree(num),Degree(denom));
    poly1 := &+[ Coefficient(num,i)*var1^i*var2^(degdeg-i) : i in [0..Degree(num)]];
    poly2 := &+[ Coefficient(denom,i)*var1^i*var2^(degdeg-i) : i in [0..Degree(denom)]];
    PrintFile(filename,Sprintf("map%oto%o := map<X%o -> X%o | [%o,%o]>;",subnum,cover,subnum,cover,poly1,poly2));
  end if;  
  if (gen eq 1) then
    var1 := Xsub.1;
    var2 := Xsub.2;
    var3 := Xsub.3;
    PrintFile(filename,Sprintf("%o := X%o.1;",var1,subnum));
    PrintFile(filename,Sprintf("%o := X%o.2;",var2,subnum));
    PrintFile(filename,Sprintf("%o := X%o.3;",var3,subnum));
    poly1 := Evaluate(s2[1],[var1,var2,var3]);
    poly2 := Evaluate(s2[3],[var1,var2,var3]);
    PrintFile(filename,Sprintf("map%oto%o := map<X%o -> X%o | [%o,%o]>;",subnum,cover,subnum,cover,poly1,poly2));
  end if;
  if (gen eq 2) then
    var1 := Xsub.1;
    var2 := Xsub.2;
    var3 := Xsub.3;
    PrintFile(filename,Sprintf("%o := X%o.1;",var1,subnum));
    PrintFile(filename,Sprintf("%o := X%o.2;",var2,subnum));
    PrintFile(filename,Sprintf("%o := X%o.3;",var3,subnum));
    defin := DefiningPolynomials(Inverse(phi));
    poly1 := Evaluate(defin[1],[var1,var2,var3]);
    poly2 := Evaluate(defin[3],[var1,var2,var3]);
    PrintFile(filename,Sprintf("map%oto%o := map<X%o -> X%o | [%o,%o]>;",subnum,cover,subnum,cover,poly1,poly2));
  end if;
  if (gen ge 3) then
    var1 := Xsub.1;
    var2 := Xsub.2;
    var3 := Xsub.3;
    PrintFile(filename,Sprintf("%o := X%o.1;",var1,subnum));
    PrintFile(filename,Sprintf("%o := X%o.2;",var2,subnum));
    PrintFile(filename,Sprintf("%o := X%o.3;",var3,subnum));
    PrintFile(filename,Sprintf("map%oto%o := map<X%o -> X%o | [%o,%o,%o]>;",subnum,cover,subnum,cover,var1,var2,var3));    
  end if;
  PrintFile(filename,Sprintf("return map%oto%o;",subnum,cover));
end if;

// Step 7c - Write the modular functions to a file if necessary

if (inftyratpoints eq true) and (subnum in covergpset) then
  fileout := "modfunc" cat IntegerToString(subnum) cat ".txt";
  printf "Writing Fourier expansion of generator to file.\n";
  System("rm " cat fileout);
  KK<z> := CyclotomicField(N);
  RRRRR<qq> := PuiseuxSeriesRing(KK);
  PrintFile(fileout,Sprintf("KK<z> := CyclotomicField(%o);\n",N));
  PrintFile(fileout,"RRRRR<qq> := PuiseuxSeriesRing(KK);");
  if (gen eq 0) then
    expdenom := ExponentDenominator(newhaupt);
    bestexpdenom := LCM([ Denominator(a/expdenom) : a in [expdenom*Valuation(newhaupt)..expdenom*AbsolutePrecision(newhaupt)-1] | Coefficient(newhaupt,a/expdenom) ne 0]);
    endprec := Floor(bestexpdenom*AbsolutePrecision(newhaupt))/bestexpdenom;
    reparsedhaupt := BigO(qq^endprec);
    for m in [bestexpdenom*Valuation(newhaupt)..bestexpdenom*endprec-1] do
      cof := KK!Eltseq(Coefficient(newhaupt,m/bestexpdenom));
      reparsedhaupt := reparsedhaupt + cof*qq^(m/bestexpdenom);
    end for;
    PrintFile(fileout,"haupt := ");
    PrintFileMagma(fileout,reparsedhaupt);
    PrintFile(fileout,";");
    PrintFile(fileout,"return haupt;");
  end if;
  if (gen eq 1) then
    expdenom := ExponentDenominator(xmodfunc);
    bestexpdenom := LCM([ Denominator(a/expdenom) : a in [expdenom*Valuation
(xmodfunc)..expdenom*AbsolutePrecision(xmodfunc)-1] | Coefficient(xmodfunc,a/expdenom) ne 0]);
    endprec := Floor(bestexpdenom*AbsolutePrecision(xmodfunc))/bestexpdenom;
    reparsedxmodfunc := BigO(qq^endprec);
    for m in [bestexpdenom*Valuation(xmodfunc)..bestexpdenom*endprec-1] do
      cof := KK!Eltseq(Coefficient(xmodfunc,m/bestexpdenom));
      reparsedxmodfunc := reparsedxmodfunc + cof*qq^(m/bestexpdenom);
    end for;
    PrintFile(fileout,"xcoord := ");
    PrintFileMagma(fileout,reparsedxmodfunc);
    PrintFile(fileout,";");
    expdenom := ExponentDenominator(ymodfunc);
    bestexpdenom := LCM([ Denominator(a/expdenom) : a in [expdenom*Valuation
(ymodfunc)..expdenom*AbsolutePrecision(ymodfunc)-1] | Coefficient(ymodfunc,a/expdenom) ne 0]);
    endprec := Floor(bestexpdenom*AbsolutePrecision(ymodfunc))/bestexpdenom;
    reparsedymodfunc := BigO(qq^endprec);
    for m in [bestexpdenom*Valuation(ymodfunc)..bestexpdenom*endprec-1] do
      cof := KK!Eltseq(Coefficient(ymodfunc,m/bestexpdenom));
      reparsedymodfunc := reparsedymodfunc + cof*qq^(m/bestexpdenom);
    end for;
    PrintFile(fileout,"ycoord := ");
    PrintFileMagma(fileout,reparsedymodfunc);
    PrintFile(fileout,";");
    PrintFile(fileout,"return xcoord, ycoord;");    
  end if;
end if;

// Step 7d - Update the status files

printf "Inftyratpoints = %o.\n",inftyratpoints;
if not fileexists("doneset.txt") then
  PrintFile("doneset.txt","dontdoset := {};");
  PrintFile("doneset.txt","computedset := {};");
  PrintFile("doneset.txt","return dontdoset, computedset;");
end if;
str := Read("doneset.txt");
dontdoset, computedset := eval str;
Include(~computedset,subnum);
if (inftyratpoints eq false) then
  queue := [subnum];
  while #queue gt 0 do
    cur := queue[1];
    Remove(~queue,1);
    if (cur ne subnum) then
      Include(~dontdoset,cur);
    end if;
    for mmmmm in [1..#newsublistmaxsub[cur]] do
      subsubsub := newsublistmaxsub[cur][mmmmm];
      if (not (subsubsub in queue)) and (not (subsubsub in dontdoset)) then
        Append(~queue,subsubsub);
      end if;
    end for;
  end while;
end if;
System("rm doneset.txt");
PrintFile("doneset.txt",Sprintf("dontdoset := %o;",dontdoset));
PrintFile("doneset.txt",Sprintf("computedset := %o;",computedset));
PrintFile("doneset.txt","return dontdoset, computedset;");
quit;