BOUND := 10000;

SetColumns(0); // No line breaks!

out := Open("attemptC.data","w");
fprintf out, 
       "ramset := SequenceToList([0:i in [1..1000]]);\nbadslopes:=SequenceToList([0:i in [1..1000]]);\n";
Flush(out);

function charpoly(A)
// Fast, non-proved, charpoly.
   if Type(A[1,1]) eq FldRatElt then
      return CharacteristicPolynomial(A : Al := "Modular", Proof:=false);
   else
      return CharacteristicPolynomial(A);
   end if;
end function;

for k:=132 to 1000 by 2 do
   t := Cputime();
   print "// +++weight",k,".";
   print "// +++Computing space of symbols";
   Mk := ModularSymbols(1,k,-1);
   M  := CuspidalSubspace(Mk);
   printf "// +++(time = %o seconds)\n", Cputime(t);
   max:= 0;
   ramset := [];
   badset := [];
   badslopes := [];
   if Dimension(M) le 1 then
      print "// +++Dimension too small";
   else
      print "// +++computing T_2";
      T := Tn(M,2);
      print "// +++computing char poly";
      f := charpoly(T);
      N := NumberField(f);
      O := EquationOrder(f);
      print "// +++computing discriminant";
      d := Discriminant(O);
      print "// +++computing primes.";
//      pset := [p:p in [1..10000] | IsPrime(p) and d mod p eq 0];
      pset := [p:p in [1..BOUND] | IsPrime(p) and d mod p^2 eq 0];
      print "// +++Looping over",pset;
      t0 := Cputime();
      for p in Reverse(pset) do
         O := pMaximalOrder(O,p);
         print "// +++",p;
      end for;
      printf "// +++Got maximalish order (%o seconds)\n", Cputime(t0);

      print "// +++Computing its disc.";
      dnew := Discriminant(O);
//      ramset := [p:p in [1..10000] | IsPrime(p) and dnew mod p eq 0];
      ramset := [p:p in [1..BOUND] | IsPrime(p) and dnew mod p^2 eq 0];

      badset:=[p : p in ramset | dnew mod p^2 eq 0];
      if #badset eq 0 then
         print "// +++No bad primes.";
      else
         printf "// +++Bad primes at weight k=%o: %o\n",k,badset;
      end if;
      print "// +++Computing slopes.";
  
      badslopes := 
        [<p,NewtonSlopes(charpoly(Tn(M,p)),p)> : p in badset];
      if #badset gt 0 then
         max := Maximum([Maximum([Denominator(Rationals()!q):q in y[2]]) : y in badslopes]);
      end if;
   end if;
   fprintf out,"badset[%o]:=%o;\n",k,badset;
   fprintf out,"badslopes[%o]:=%o;\n",k,badslopes;
   Flush(out);
   printf "badset[%o]:=%o;\n",k,badset;
   printf "badslopes[%o]:=%o;\n",k,badslopes;
   if max gt 1 then
      print  "// +++max:=",max;
   end if;
   if max gt 2 then
      print "// +++HOORAY!";
   end if;

   print "// +++ time = ",Cputime(t);

   DisownChildren(Mk);
end for;

/*
***

I know, I know, it's a long program. but it was written by me so should easily
be understandable because I only know a few commands. What it does is loops
through k, computing modular symbols M of level 1 and weight k, and then
computes the number field N generated by a newform. The question I am
desperate to answer is whether I can find a prime p for which the
p-adic valuation of a_p has denominator >=3. I am working on the
assumption that if p>10000 then I'm unlikely to be able to work out
T_p anyway, so I'm restricting to small p, hoping that the first
example, if there is one, will be small. So I compute O, the equation
order for N (N has been defined using the char poly f of T_2) and
then the small primes dividing disc(f). These are the primes that potentially
ramify in N. To check they don't I try to bump up O so it's maximal at
each such prime (David: I should loop through Reverse(pset) instead of pset,
for some reason; it speeds the whole process up---do you know why? See
earlier email). Because at the end of the day I'm looking for something
with valuation with denominator >=3, I need N to have a prime ideal
with "e>=3" i.e. a prime \p of N such that (\p)^3 divides p, its
res char. This I think implies that p will divide the discriminant
at least twice. ramset is the small ramified primes, badset is the ones
dividing disc(O) at least twice (once O has been bumped up).
Because Decomposition(O,p) is broken in the version I have (David,
I told you about this: if O is maximal at p it should work, right?) , I then
have to compute T_p. So I check. I print the results to a file.

None of the above is relevant, but it's what I happen to be doing
in my spare time. What is relevant is that it seems to me that
when the loop finishes, I need to remember *none* of the variables
(except for the output file variable);
the only things I'm interested in are ramset and badset (and badslopes, out
of interest, given that I calculated them), but I've fprintf'ed them
out so I can even kill them. So I simply cannot understand
why the amount of memory that this loop takes just goes up and up and up.
It doesn't take long to hit 10M, according to the top command,
and if I leave it on overnight it gets to 60 meg. 

I pestered William to find out how Disownchildren works, and I think
I'm using it correctly. I know I've sent both of you a lot of emails
recently, for which apologies (although David: I hope at least a few
have been helpful!!)---but I still simply don't understand why the
amount of memory used by this loop doesn't periodically drop! I have
found this really funny trick: run "top" and then put "s0(return)"
and top just goes mental! But it seems to indicate that the amount
of memory used by magma just goes up and up and up. I really don't
want this to happen because swapping shouldn't be inevitable in
this loop. Not at weight 200, anyway...

Kevin

PS this email finally isolates the thing I don't understand, so
rest assured I won't pester you much more with this. Ooh, my laptop
screensaver has just kicked in and it's those amazing 3d sharks.
I love them.


*/
