function d(k)
   return DimensionCuspFormsGamma0(1,k);
end function;

intrinsic BuzzardSlopes(kmax::RngIntElt) -> SeqEnum
{}
   return BuzzardSlopes(kmax,[]);
end intrinsic;

intrinsic BuzzardSlopes(kmax::RngIntElt, slopes::SeqEnum) -> SeqEnum
{The slopes of U_2 on weight-k forms for SL_2(Z), for k=2,4,6,...,kmax.  
The first argument "slopes" is the result of a previous invocation,
and is required only for reasons of efficiency.}

   if #slopes eq 0 then
      slopes := [[]];
   end if;   
   kstart := Max([4] cat [i : i in [1..#slopes] | slopes[i] ne []]);
   slopes := slopes cat [[] : i in [1..kmax - #slopes]];
   
   for k := kstart to kmax by 2 do
      n := 0;
      while 2^(n+1) lt k do
         n +:= 1;
      end while;
      gamma := k - 2^n;
      beta  := gamma - 1;
      alpha := (beta-1)/2;
      w0    := slopes[gamma];
      v1    := slopes[k-2^(n-1)];
      v2    := [v1[i+d(gamma)] : i in [1..#v1 - d(gamma)]];
      v     := [];
      xx    := [];
      for i in [1..#v2] do
         if v2[i] lt alpha then
            Append(~v,v2[i]);
         end if;
         if v2[i] eq alpha then
            Append(~xx,alpha);
         end if;
      end for;
      w0 := w0 cat [v[i]+1 : i in [1..#v]]
            cat xx 
            cat [2*alpha-1-v[#v+1-i] : i in [1..#v]]
            cat [beta - slopes[gamma][#slopes[gamma]+1-i] : 
                  i in [1..#slopes[gamma]] ];
      if gamma gt 2 then
         Append(~w0, beta);
      end if;

      w0 := w0 cat [slopes[k-2*beta][i]+beta :
                        i in [1..#slopes[k-2*beta]]];

      slopes[k] := [w0[i] : i in [1..d(k)]];
   end for;

   return slopes;
end intrinsic;

intrinsic Level2Slopes(kmax::RngIntElt) -> SeqEnum
{The slopes of U_2 on weight-k forms for Gamma_0(2), for k=2,4,6,...,kmax.}
   return Level2Slopes(kmax,[]);
end intrinsic;

function ExtendToLevel2(k, level1slopes)
   // The slopes of U_2 on weight-k forms for Gamma_0(2), for k=2,4,6,...,kmax.
   if IsOdd(k) then
      return level1slopes;
   end if;
   for i in [1..#level1slopes] do
      Append(~level1slopes, k-1-level1slopes[i]);
   end for;
   for i in [1..DimensionNewCuspFormsGamma0(2, k)] do
      Append(~level1slopes, (k div 2) - 1);
   end for;
   return Sort(level1slopes);   
end function;

intrinsic Level2Slopes(kmax::RngIntElt, level1slopes::SeqEnum) -> SeqEnum, SeqEnum
{The slopes of U_2 on weight-k forms for Gamma_0(2), and the slopes for level 1,
for k=2,4,6,...,kmax.   The last argument, slopes, is the slopes of level1, which
is the second argument returned by this intrinsic.}
   /*
   1. Call Buzzard's level 1 slopes routine to compute the level 1 
      slopes for weight up to kmax. 
   2. Add level 2 slopes onto the list of level 1 slopes computed 
      in step 2.
   */
   level1slopes := BuzzardSlopes(kmax, level1slopes);
   return [ExtendToLevel2(k, level1slopes[k]) : k in [1..#level1slopes]], 
          level1slopes;
end intrinsic;

intrinsic MazurS(alpha::RngIntElt, 
                     N::RngIntElt, 
                     p::RngIntElt,
                slopes::SeqEnum) -> FldRatElt
{}
   return MazurS(1.0*alpha,N,p,slopes);
end intrinsic;


intrinsic MazurS(alpha::FldPrElt, 
                     N::RngIntElt, 
                     p::RngIntElt,
                slopes::SeqEnum) -> FldPrElt
{Let d(k) be the number of newforms
on SL_2(Z) of weight k, and s(k,alpha) be the number of these newforms
having supersingularity <= alpha.  The supersingularity of a newform
of weight k is the slope of the newform divided by p+1.
MazurS is the rational number 1/N * sum_\{k<=N\} s(k,alpha)/d(k).
}
   require N ge 1 : "Argument 2 must be at least 1.";
   require IsPrime(p) : "Argument 3 must be prime.";

   sum := 0.0;
   cnt := 0;
   for k in [i : i in [2..N] | IsEven(i) and #slopes[i] gt 0] do
      sum +:= #[s : s in slopes[k] | s le (k-1)*alpha]*1.0 / #slopes[k];
      cnt +:= 1;
   end for;
   return 1/cnt * sum;
end intrinsic;

