/***********************************************************
 j0n_invariants.m -- compute invariants of factors of J_0(N)

    For each level N, compute lots of information about 
    the simple factors of J_0(N).

    William Stein, 28 January 2001.

 ***********************************************************/

forward OddPart,  
        PrimeSeq;

J0N := recformat<
   levels,
   new_dimensions,
   dimensions,
   intersection_graph,
   ap_traces,
   hecke_fields,
   atkin_lehners,
   component_group_orders,
   tamagawa_numbers,
   torsion_upper_bounds,
   torsion_lower_bounds,
   l_ratios,
   analytic_sha_upper_bounds,
   analytic_sha_lower_bounds,
   eigenvalues,
   q_expansions
>;


intrinsic J0(N::RngIntElt : Sign := +1, bound := 38) -> Rec
{Lots of data about the factors of J_0(N).}
   M := ModularSymbols(N,2,Sign);
   S := CuspidalSubspace(M);
   D := NewformDecomposition(S : Proof := false);
   D := SortDecomposition(D);
   R<x> := PolynomialRing(Rationals());

   intersection_graph := 
                      [OddPart(#IntersectionGroup(A,B) cmpeq Infinity() 
                            select 0 else #IntersectionGroup(A,B), Sign) : 
                                                            A in D, B in D];

   ap_traces := [[Trace(DualHeckeOperator(A,p)) : p in PrimeSeq(2,bound)]
   		   : A in D | IsNew(A)];

   atkin_lehners := [[AtkinLehner(A,p[1])[1,1] : p in Factorization(N)]
     		        : A in D | IsNew(A)];

   q_expansions := [* *];
   eigenvalues := [* *];
   for A in D do
      if IsNew(A) then
         f := qEigenform(A,bound);
         Append(~q_expansions, f) ;
         Append(~eigenvalues, [Coefficient(f,p) : p in PrimeSeq(2,bound)]);
      end if;
   end for;
   hecke_fields := [];
   for f in q_expansions do
      if BaseRing(Parent(f)) cmpeq RationalField() then
         Append(~hecke_fields, x-1);
      else
         Append(~hecke_fields, R!Modulus(BaseRing(Parent(f))));
      end if;
   end for;

   square_free := [p[1] : p in Factorization(N) | p[2] eq 1];
   component_group_orders := [[OddPart(ComponentGroupOrder(A,p), Sign) : p in square_free]
                                 : A in D | IsNew(A)];
   tamagawa_numbers := [[OddPart(TamagawaNumber(A,p), Sign) : p in square_free]
                                 : A in D | IsNew(A)];
   torsion_upper_bounds := [OddPart(TorsionBound(A,bound), Sign) : A in D |IsNew(A)];
   e := M!<1,[Cusps()|0,Infinity()]>;
   torsion_lower_bounds := [OddPart(#SubgroupOfTorus(A,e), Sign) : A in D |IsNew(A)];

   l_ratios := [OddPart(LRatio(A,1), Sign) : A in D | IsNew(A)];
 
   analytic_sha_lower_bounds := 
          [(l_ratios[i]*torsion_lower_bounds[i]^2)/&*tamagawa_numbers[i] : 
              i in [1..#l_ratios]];
   analytic_sha_upper_bounds := 
          [(l_ratios[i]*torsion_upper_bounds[i]^2)/&*tamagawa_numbers[i] : 
              i in [1..#l_ratios]];

   ans := rec<J0N | 
      levels := [Level(AssociatedNewSpace(A)) : A in D],
      dimensions := [Dimension(A) : A in D],
      new_dimensions := [Dimension(AssociatedNewSpace(A)) : A in D],
      intersection_graph := intersection_graph,
      ap_traces := ap_traces,
      hecke_fields := hecke_fields,
      atkin_lehners := atkin_lehners,
      component_group_orders := component_group_orders,
      tamagawa_numbers := tamagawa_numbers,
      torsion_lower_bounds := torsion_lower_bounds,
      torsion_upper_bounds := torsion_upper_bounds,
      l_ratios := l_ratios,
      analytic_sha_upper_bounds := analytic_sha_upper_bounds,
      analytic_sha_lower_bounds := analytic_sha_lower_bounds,
      q_expansions := q_expansions,
      eigenvalues := eigenvalues
   >;

   DisownChildren(M);

   return ans;
end intrinsic;


function OddPart(n, sign)
   if sign eq 0 then
      return n;
   end if;
   return n eq 0 select 0 else Parent(n)!(n / 2^Valuation(n,2));
end function;

function PrimeSeq(n,m)
   return [p : p in [n..m] |IsPrime(p)];
end function;


function Go(N, file)
   return Sprintf("echo \"SetVerbose(\\\\\"ModularForms\\\\\",2); Attach(\\\\\"j0n_invariants.m\\\\\"); t := Cputime(); a := J0(%o); fprintf Open(\\\\\"%o\\\\\", \\\\\"a\\\\\"), \\\\\"\\\\n\\\\nJ[\\%o] := \\%o ;  // time = \\%o seconds\\\\n\\\\n\\\\n\\\\\", %o, a, Cputime(t); quit;\" | magma\n", N, file, N);
end function;

intrinsic MakeGoFile(Nstart, Nstop)
{}
   fname := "go_" cat IntegerToString(Nstart) cat "-" cat IntegerToString(Nstop);
   out   := "j0_" cat IntegerToString(Nstart) cat "-" cat IntegerToString(Nstop) cat ".out";
   file := Open(fname, "w");
   for n in [m : m in [Nstart..Nstop] |IsSquareFree(m)] do 
      fprintf file, "%o\n", Go(n,out); 
   end for;
end intrinsic;

intrinsic IsSquareFree(N) -> BoolElt
{}
   return &*[Integers()|p[2] : p in Factorization(N)] eq 1;
end intrinsic;

// for n in [m : m in [1..500] |IsSquareFree(m)] do Go(n,"1-500.out"); end for;
