/***********************************************************************
NOTE: 
  * I should write a p-adic characters package. -- WAS
  * His function name "Lseries" is bad.   It should be
    pAdicLSeriesApproximation
 ***********************************************************************/

/*
Date: Mon, 29 Jan 2001 01:26:54 -0500 (EST)
From: Robert Pollack <rpollack@math.harvard.edu>
To: "William A. Stein" <was@math.harvard.edu>
Subject: Lseries code

> Hey, send me some magma code that involves p-adic L-functions!

Okay...but the code is really awful.  I believe it works though.

I will send you two more e-mails.  One containing a write-up in
TeX of computing the p-adic L-series and the other will be my
code.  The first will probably be more helpful than the later.

As the code is not documented, let me at least lead you thru 
what it is supposed to do.  The main procedure is called Lseries.
It takes in five parameters - let's call them V,p,a_p,n,D.

The V should be the space of modular symbols associated to E
The p is the prime
The a_p is the a_p value of E
The n represents the level of accuracy.
The D represents twisting by a quadratic character of conductor D.
         (you can always just take D=1)

Then Lseries(V,p,a_p,n,D) returns in the ordinary case a polynomial
of degree p^(n-1)-1 approximating the L-series.  In the supersingular
case, it returns two polynomials say G and H so that G + H * alpha
approximates the L-series where alpha is a root of x^2 - a_p x + p.

The other main procedure is iwasawa_invariants_w_curve.  This
takes 6 arguments, say N,V,p,a_p,n,D.

Here N is the conductor of the curve
     V is the space of modular symbols of the curve
     p is some prime
     a_p is the a_p value of the curve
     n is the level of accuracy
     D represents twisting by a quadratic character of conductor D

This procedure is supposed to return the iwasawa invariants of 
the L-series.  

Precisely, in the ordinary case it returns a 6-tuple of data,
namely the rank of the curve, 
       the p-adic valuation of the leading non-zero term of the L-series,
       the mu invariant of the L-series
       the lambda invariant of the L-series
       the newton polygon of the L-series.
                 -this is returned as [ [m_1,s_1] , [m_2,s_2] , ... ]
                 where m_1 is the number of roots of slope s_1...;
                 a slope of 10000 represents the root 0.
       a boolean saying whether the above data is correct.
         (note if mu is positive then a "true" only means that the data
         is as accurate as it can be using a polynomial of degree
         p^(n-1)-1; that is it could have a another root of smaller slope)


In the supersingular case, the p-adic L-series is not an iwasawa
function so the above data couldn't make any sense.  Instead it
returns an 11-tuple of data consisting of the rank and then the
same data for the g and h functions i construct in my thesis.  i should
have a rough draft of my thesis ready soon if you are interested
in how they are defined.  basically, they are the "real" and "imaginary"
parts of the L-series with their trivial zeroes removed.  They are
in fact iwasawa functions (actually integral power series).

I've been back and forth between NY and Boston for the last
few weeks.  Now I should be back in Boston full time for the
next month.

Okay...let me try and send out the code.

Robert
*********************************************************************/

R<x>:=PolynomialRing(Rationals());
R<T>:=PolynomialRing(Rationals());

//This function cuts out the subspace of H_1(X_0(N),Z) where
//T_2 acts as a_2, T_3 acts as a_3, T_5 acts as a_5 and T_7 acts as a_7
//Input:  	N - conductor
//		a_2 - eigenvalue of T_2
//		a_3 - eigenvalue of T_3
//		a_5 - eigenvalue of T_5
//		a_7 - eigenvalue of T_7
//
//Output:	space of modular symbols of level N and specified eigenvalues
function form_curve(N,a_2,a_3,a_5,a_7)

   R<x>:=PolynomialRing(Rationals());
   M:=ModularSymbols(N,2,+1);
   I:=[<2,x-a_2>,<3,x-a_3>,<5,x-a_5>,<7,x-a_7>]; 
   V:=Kernel(I,M);
   ModularSymbolEven(V,1/2);

   return(V);
end function;

function table_of_modular_symbols(V,p,n,D) 
   local q,table,temp,k,scale,q1,mid,temp2,e;

   if (D mod p eq 0) then
      print "Twisting by something not relatively prime to p!";
   end if;

   if (KroneckerSymbol(D,-1) eq -1) then
      print "Twisting by an odd character!";
   end if;

   q1 := p^(n-1);
   q := q1*p;

   k:=1;
   while (ModularSymbolEven(V,k/p^10)[1] eq 0) do
      k := k+1;
   end while;

   scale := Denominator(ModularSymbolEven(V,k/q)[1]);
   table := [];

   for a in [0..(q-1) div 2] do
      if (Valuation(a,p) le 2) or (a eq 0) then
         temp2 :=0;
         for b in [1..D] do
            e := KroneckerSymbol(D,b);
            if (e ne 0) then
               temp := ModularSymbolEven(V,(D*a-q*b)/(D*q))[1] * scale;
               if Denominator(temp) ne 1 then   
                  print "PROBLEMS!  There are still denominators in ",
                        "the modsyms.";
               end if;   
               temp2 := temp2 + KroneckerSymbol(D,b)*temp;
            end if;
         end for;
         table := table cat [Numerator(temp2)];
      else
         table := table cat [0];
      end if;
   end for;

   //USING the fact that [r/s]^+ equals [-r/s]^+
   mid := ((q-1) div 2) + 1;
   for a in [mid..q-1] do
      table := table cat [table[q-a+1]];
   end for;
   table := table cat [scale];

   return(table);
end function;

//Multiplies a and b interpreted as a[1] + sqrt(-p)*a[2].
function multiply(a,b,p)
   return([a[1]*b[1] - p*a[2]*b[2],a[1]*b[2]+a[2]*b[1]]);
end function;

//Here answer[a] = [alpha^-a]
//n+1 is the maxl power and m is the accuracy
function powers_of_alpha(p,a_p,n,m)
   local temp,answer,t,q;

   if a_p mod p ne 0 then
      temp := a_p mod p;

      for k in [2..m] do
         t := -(temp^2-a_p*temp+p) div p^(k-1) mod p *
            InverseMod(2*temp - a_p,p);
         temp := temp + t*p^(k-1);
      end for;

      q := p^m;
      temp := InverseMod(temp,q);
      answer := [[temp]];

      for k in [2..n+1] do
         Append(~answer,[answer[k-1][1]*answer[1][1] mod q]);
      end for;

      return(answer);
   else
      if a_p eq 0 then
         temp := [0,-1/p];
      else
         if a_p eq -3 then
            temp := [-1/2,-1/6];
         else
            temp := [1/2,-1/6];
         end if;
      end if;
      answer := [temp];
      for k in [2..n+1] do
         Append(~answer,multiply(answer[k-1],answer[1],p));
      end for;
      return(answer);
   end if;
end function;

//Here answer = {a} mod p^m
function one_teich(a,p,m)
   local answer,t;

   answer := a mod p;
   for k in [2..m] do
      t := 0;
      while (answer + t*p^(k-1))^(p-1) mod p^k ne 1 do
         t := t+1;
      end while;
      answer := answer + t*p^(k-1); 
   end for;

   return(answer);
end function;
         
//Here answer[a] = {a} (mod p^m) 
function list_of_teich(p,m)
   local answer;

   answer := [];
   for a in [1..p-1] do
      Append(~answer,one_teich(a,p,m));
   end for;
   return(answer);
end function;

//Here answer[a] = gamma^(a-1)
//SO THERE IS A SHIFT
//n is the maxl power and m is the accuracy
function powers_of_gamma(gamma,p,n)
   local answer;

   answer := [1];
   for k in [2..p^(n-1)] do
      Append(~answer,answer[k-1]*gamma mod (p^n));
   end for;

   return answer;
end function;


function lift_polynomial(f)
   return PolynomialRing(Rationals()) ! 
                (PolynomialRing(Integers()) ! f);
end function;

//Here answer[a]=(1+T)^(a-1)
//SO THERE IS A SHIFT
//n is the maxl power and m is the accuracy
function powers_of_oneplusT(p,n,m)
   local R,S,answer,q;

   q := p^m;
   R := ResidueClassRing(q);
   S<T> := PolynomialRing(R);
   answer := [1];
   for k in [2..p^(n-1)] do
      answer := answer cat [answer[k-1]*(1+T)] ;
   end for;

   return(answer);
end function;

function alpha(n,ALPHA)
     return(ALPHA[n]);
end function;

function modsym(a,MODSYM)
   return(MODSYM[a+1]);
end function;

function oneplusT(j,POWER)
   return(POWER[j+1]);
end function;

function teich(a,TEICH)
   return(TEICH[a]);
end function;

function gamma(j,GAMMA)
   return(GAMMA[j+1]);
end function;

function scale(list,c)
   local k,answer;

   answer:=[];
   for k in [1..#list] do
      answer := answer cat [list[k]*c];
   end for;
   return(answer);
end function;

function add(list1,list2)
   local k,answer;
   
   if #list1 ne #list2 then
      print "Error in add";
   end if;

   answer:=[];
   for k in [1..#list1] do
      answer := answer cat [list1[k]+list2[k]];
    end for;

   return(answer);
end function;

function subtract(list1,list2)
   local k,answer;
   
   if #list1 ne #list2 then
      print "Error in subtract";
   end if;

   answer:=[];
   for k in [1..#list1] do
      answer := answer cat [list1[k]-list2[k]];
    end for;

   return(answer);
end function;

function measure(a,p,n,ALPHA,MODSYM)
   local q,q1;

   q := p^n; 
   q1 := q div p;
   return( subtract(scale(alpha(n,ALPHA),modsym(a,MODSYM)),
           scale(alpha(n+1,ALPHA),modsym((a mod q1) * p ,MODSYM))));

end function;

//This procedure returns a polynomial approximation to the p-adic
//L-series of the elliptic curve corresponding to V twisted by a quadratic
//character of discriminant D.  The degree of the polynomial is p^n-1. 
//If a_p is prime to p (ordinary) then the answer is simply a polynomial.
//Otherwise (in the supersingular case) the answer is in the form [G,H]
//where the L-series is approximated by G + sqrt(-p) * H.
//
//Input:   V - space of modular symbols
//      p - prime
//      a_p - eigenvalue of T_p acting on V
//      n - relates to degree of approximation
//      D - value to twist by
function Lseries(V,p,a_p,n,D)
   local answer,q;

   ALPHA := powers_of_alpha(p,a_p,n,n+1);
   MODSYM := table_of_modular_symbols(V,p,n,D);
   POWER := powers_of_oneplusT(p,n,n+1);
   TEICH := list_of_teich(p,n+1);
   GAMMA := powers_of_gamma(1+p,p,n+1); 

   q1 := p^(n-1);
   q := q1*p;
   if a_p mod p ne 0 then
      answer := [0];
   else
      answer := [0,0];
   end if;

   for a in [1..p-1] do
      for j in [0..q1-1] do
         answer := add(answer, 
                   scale(measure(teich(a,TEICH)*gamma(j,GAMMA) mod 
                                q,p,n,ALPHA,MODSYM),
                        lift_polynomial(oneplusT(j,POWER))));
      end for;
   end for;

   answer:=scale(answer,1/D);
   answer:=scale(answer,1/MODSYM[#MODSYM]);

   return(answer);
end function;

function highest_power(n,p)
   local k,temp;

   k:=0;
   temp:=1;
   while (n ge temp) do
      temp := temp*p;
      k := k+1;
   end while;
   return(k-1);
end function;

function valuation_of_coefficients(f,p,n)
   local answer,j,temp,okay,bound,list_bad;

   answer := [];
   okay:=true;
   list_bad:=[];
   for j in [1..Degree(f)+1] do
      if Coefficient(f,j-1) ne 0 then
         answer := answer cat [Valuation(Coefficient(f,j-1),p)];
      else 
         answer := answer cat [10000];
      end if;
      
      bound := n-1-highest_power(j-1,p);
      if (Valuation(Coefficient(f,j-1),p) ge bound) and (j ne 1) then
         list_bad := list_bad cat [j];
      end if;
   end for;

   if f eq 0 then
      answer := [];
      answer := answer cat [10000] cat [10000] cat [10000];
   end if;
   
   return(<answer,list_bad>);
end function;
   
function newton_polygon(C)
   local k,here,next,min_slope,answer,okay,prev_slope;

   answer:= [];
   okay:=true;
   here:=1;
   next:=1;

   while here lt #C do
      min_slope:=100000;
      for k in [here+1..#C] do
         if (C[k]-C[here])*(1.0)/(k-here) le 1.0*min_slope then
            min_slope := (C[k]-C[here])/(k-here);
            next := k;
         end if;
      end for;
      answer := answer cat [[next-here,-min_slope]];
      here := next;
   end while;   
   return(answer);
end function;

function find_last_nontriv(NP,p,sign)
   local slope,place;

   place := 0;
   if sign mod 2 eq 0 then
      slope := 1/(p*(p-1));
   else
      slope := 1/(p-1);
   end if;

   for k in [1..#NP] do
      if (NP[k][2] gt 0) then
         if (NP[k][2] ne slope) or (NP[k][1] ne 1/slope) then
            place:=k;
         else
            slope:=slope/(p^2);
         end if;
      end if;
   end for;

   return(place);
end function;


function refine_newton_polygon(NP,p,a_p,sign)
   local answer,k,last;
   
   answer := [];
   if a_p mod p ne 0 then
      for k in [1..#NP] do
         if NP[k][2] gt 0 then
            answer := answer cat [NP[k]];
         end if;
      end for;
   else
      last := find_last_nontriv(NP,p,sign);
      for k in [1..last] do
         answer := answer cat [NP[k]];
      end for;
   end if;
            
   return(answer);
end function;

function minimize_ord(C,p,n) 
   local k,start;

   start:=1;
   if C[1] eq 10000 then
      start := 2;
   end if;

   for k in [start..#C] do
      if C[k] ge n-1-highest_power(k-1,p) then
         C[k] := n-1-highest_power(k-1,p);
      end if;
   end for;
   return(C);
end function;

function in_list(list,num)
   local k,flag;

   flag:=false;
   for k in [1..#list] do
      if list[k] eq num then
         flag := true;
      end if;
   end for;
   
   return(flag);
end function;

function compute_rank(C,N,r,D)
   local rank;

   if C[1] eq 10000 then
      rank := 1;
   else
      rank :=0;
    end if;

   if (r eq 0) and (rank eq 1) and (KroneckerSymbol(D,N) eq 1) then  
          //it has even rank >= 2.
      rank:=2;
      C[2]:=10000;
   end if;

   if (r eq 1) and (rank eq 1) and (KroneckerSymbol(D,N) eq -1) then  
          //it has even rank >= 2.
      rank:=2;
      C[2]:=10000;
   end if;

   return(<C,rank>);
end function;


function check_NP(NP,C,list_bad,rank,p,a_p,n,sign)
   local temp,val,k,good;

   //checks if break points are accurate.
   
   good:=true;
   temp := refine_newton_polygon(NP,p,a_p,sign);
   val := 1;
   for k in [1..#temp] do
      val := val + temp[k][1];
      if (in_list(list_bad,val)) then
         good := false;
         print "Failed at a breakpoint,",val;
      end if;
   end for;

   //checks if lowering the valuation of any unknown coefficient to its
   //minimal possible value affects the newton polygon.

   C_min := minimize_ord(C,p,n); //lowers all unknown coefficients 
                                 // to lowest possible valuation
   
   if rank eq 2 then
         C_min[2] := 10000;
   end if;

   NP_min := newton_polygon(C_min);
   if refine_newton_polygon(NP_min,p,a_p,sign) 
          ne refine_newton_polygon(NP,p,a_p,sign) then
      good := false;
      print "Failed while lowering.";   
   end if;

   return(good);
end function;

function lambda(NP)
   local k,lam;

   k:=1;
   lam := 0;
   while (k le #NP) and (NP[k][2] gt 0) do
        lam := lam + NP[k][1];
      k := k+1;
   end while;

   return(lam);
end function;

function mu(NP,C,rank)
   local k,val,m,start;

   k:=1;
   val := 0;
   if (rank ne 2) then
      while (k le #NP) and (NP[k][2] gt 0) do
         if NP[k][2] lt 1000 then
            val := val + NP[k][1]*NP[k][2];
         end if;
         k := k+1;
      end while;
   else
      if (#NP gt 0) and (NP[1][1] eq 1) then 
         start := 3;
      else
         if #NP eq 0 then
            start := 1;
         else
            start := 2;
         end if;
      end if;
      k:=start;
      while (k le #NP) and (NP[k][2] gt 0) do
         if NP[k][2] lt 1000 then
            val := val + NP[k][1]*NP[k][2];
         end if;
         k := k+1;
      end while;
   end if;

   if rank eq 0 then
      m := C[1]-val;
   else
      if rank eq 1 then
         m := C[2] - val;
      else
         m := C[3] - val;
      end if;
   end if;

   return(m);

end function;

function remove_two(list)
   local answer,k;

   answer := [];

   for k in [1..#list] do
      if list[k] ne 2 then
         answer := answer cat [list[k]];
      end if;
   end for;

   return(answer);

end function;

function leading_term(C,rank)

   if rank eq 0 then
      return(C[1]);
   else
      if rank eq 1 then
         return(C[2]);
      else
         if rank eq 2 then
            return(C[3]);
         end if;
      end if;
   end if;
end function;

function compute_untwisted_rank(V,p,a_p)
   local L,C;

   L:=Lseries(V,p,a_p,2,1);
   L:=L[1];
   C := valuation_of_coefficients(L,p,1);
   if C[1][1] ne 10000 then
      return(0);
   else 
      return(1);
   end if;
end function;

function remove_unit_zeroes(NP)
   local k,answer;

   answer := [];
   for k in [1..#NP] do
      if NP[k][2] gt 0 then
         answer := answer cat [NP[k]];
      end if;
   end for;
   return(answer);
end function;

//Has the space of modular symbols already pre-computed
function iwasawa_invariants_w_curve_ord(N,V,p,a_p,n,D)
   local C,NP,k,lam,m,rank,j,list_bad,list_breapoints,
         C_min,NP_min,val,good,m2,lam2,temp,leading;

   L:=Lseries(V,p,a_p,n,D)[1];
   C := valuation_of_coefficients(L,p,n);
   list_bad:=C[2];
   C:=C[1];
      
   //FIRST COMPUTE RANK
   
   rank := compute_untwisted_rank(V,p,a_p);
   C := compute_rank(C,N,rank,D);

   rank := C[2];  //changes valuation of a_1 to 10000 if guesses rank 2.
   C := C[1];   

   if rank eq 2 then
      list_bad := remove_two(list_bad);
   end if;

   NP := newton_polygon(C);
   good := check_NP(NP,C,list_bad,rank,p,a_p,n,0);  //0 represents nothing.
   NP := remove_unit_zeroes(NP);
   leading := leading_term(C,rank);
   
   if (leading eq 0) and (in_list(list_bad,rank+1) eq false) then
        good:=true;
   end if;

   lam:=lambda(NP);
   m := mu(NP,C,rank);   

   return(<rank,leading,m,lam,NP,good>);
end function;

function scale_by_p(C)
   local k;

   for k in [1..#C] do
      if C[k] ne 10000 then
         C[k]:=C[k]+1;
      end if;
   end for;

   return(C);
end function;

   
function remove_forced_zeroes(NP,p,init_slope)
   local k,slope,temp;

   temp := [];
   k := 1;   
   slope:=init_slope;
   while (k le #NP) and (NP[k][2] gt 0)  do
      if NP[k][2] ne slope then
         temp := temp cat [NP[k]];
      else
         if NP[k][1] - 1/slope ne 0 then
            temp := temp cat [[NP[k][1] - 1/slope,slope]];
         end if;
          slope := slope/(p^2);
      end if;
      k := k+1;
   end while;

   return(temp);
end function;

//Has the space of modular symbols already pre-computed
function iwasawa_invariants_w_curve_ss(N,V,p,a_p,n,D)
   local L,G,H,k,CG,CH,mg,mh,lamg,lamh,NPG,NPH,rank,rankg,rankh,
      list_badg,list_badh,denomg,denomh,goodg,goodh,
      list_breakpointsg,list_breakpointsh,
      leadingg,leadingh;

   L := Lseries(V,p,a_p,n,D);
   G := L[1];
   H := L[2];

   //   denomg and denomh represent the number of p's 
   //   in the denominator of the measure

   if n mod 2 eq 0 then
      denomg:=n div 2;
      denomh:=(n+2) div 2;
   else
      denomg:=(n+1) div 2;
      denomh:=(n+1) div 2;
   end if;   
      
   CG := valuation_of_coefficients(G,p,n-denomg);
   CH := valuation_of_coefficients(H,p,n-denomh);

   list_badg := CG[2];
   list_badh := CH[2];

   CG:=CG[1];
   CH:=CH[1];

   //Computing the rank

   rank:=compute_untwisted_rank(V,p,a_p);

   CG := compute_rank(CG,N,rank,D); //changes a_1 if guesses rank 2
   rankg := CG[2];
   CG := CG[1];

   CH := compute_rank(CH,N,rank,D);
   rankh := CH[2];
   CH := CH[1];

   if rankg ne rankh then
      print "WEIRD..ranks different in G and H";
      PrintFile("notes","ranks different in G and H");
      PrintFile("notes",D);
   end if;

   rank := rankg;
   if rank eq 2 then
      list_badg := remove_two(list_badg);
   end if;

   if rank eq 2 then
      list_badh := remove_two(list_badh);
   end if;

   NPG := newton_polygon(CG);
   NPH := newton_polygon(CH);

   goodg := check_NP(NPG,CG,list_badg,rank,p,a_p,n,0); //0 represents g
   goodh := check_NP(NPH,CH,list_badh,rank,p,a_p,n,1); //1 represents h

   leadingg := leading_term(CG,rank);
   leadingh := leading_term(CH,rank);

   if (a_p eq 0) then
      if (leadingg eq -1) and (in_list(list_badg,rank+1) eq false) then
          goodg:=true;
      end if;

      if (leadingh eq -1) and (in_list(list_badh,rank+1) eq false) then
          goodh:=true;
      end if;
   end if;


   if a_p eq 0 then
      NPG := remove_forced_zeroes(NPG,p,1/((p-1)*p));
      NPH := remove_forced_zeroes(NPH,p,1/(p-1));
   else
      if a_p eq -3 then
         NPG := remove_forced_zeroes(NPG,p,1/((p-1)*p));
         NPH := remove_forced_zeroes(NPH,p,1/(p-1)*p^2);
      end if;
   end if;

   lamg := lambda(NPG);
   lamh := lambda(NPH);   

   mg := mu(NPG,CG,rank)+1;
   mh := mu(NPH,CH,rank)+1;

   return(<rank,leadingg,mg,lamg,NPG,goodg,leadingh,mh,lamh,NPH,goodh>);
end function;

// This procedure returns many invariants of the curve corresponding
// to V twisted by D.  
//
// In the ordinary case the data returned is:
//   <rank,valuation of leading coefficient,mu invariant,lambda invariant,
//    newton polygon of L-series, boolean saying whether the above data
//    are accurate>
//
// The newton polygon is returned as [n_1,s_1],[n_2,s_2],... where there are
// n_1 roots of slope s_1, etc.
//
// In the supersingular case the data returned is similar.  The first entry is
// the rank, then the next five are the data for my "g" iwasawa function and
// then the next five for my "h" iwasawa function.

function iwasawa_invariants_w_curve(N,V,p,a_p,n,D)
   if a_p mod p ne 0 then
      return(iwasawa_invariants_w_curve_ord(N,V,p,a_p,n,D));
   else
      return(iwasawa_invariants_w_curve_ss(N,V,p,a_p,n,D));
   end if;
end function;