attachment:lseries.m of days22/wuthrich

   1 /***********************************************************************
   2 NOTE: 
   3   * I should write a p-adic characters package. -- WAS
   4   * His function name "Lseries" is bad.   It should be
   5     pAdicLSeriesApproximation
   6  ***********************************************************************/
   7 
   8 /*
   9 Date: Mon, 29 Jan 2001 01:26:54 -0500 (EST)
  10 From: Robert Pollack <[email protected]>
  11 To: "William A. Stein" <[email protected]>
  12 Subject: Lseries code
  13 
  14 > Hey, send me some magma code that involves p-adic L-functions!
  15 
  16 Okay...but the code is really awful.  I believe it works though.
  17 
  18 I will send you two more e-mails.  One containing a write-up in
  19 TeX of computing the p-adic L-series and the other will be my
  20 code.  The first will probably be more helpful than the later.
  21 
  22 As the code is not documented, let me at least lead you thru 
  23 what it is supposed to do.  The main procedure is called Lseries.
  24 It takes in five parameters - let's call them V,p,a_p,n,D.
  25 
  26 The V should be the space of modular symbols associated to E
  27 The p is the prime
  28 The a_p is the a_p value of E
  29 The n represents the level of accuracy.
  30 The D represents twisting by a quadratic character of conductor D.
  31          (you can always just take D=1)
  32 
  33 Then Lseries(V,p,a_p,n,D) returns in the ordinary case a polynomial
  34 of degree p^(n-1)-1 approximating the L-series.  In the supersingular
  35 case, it returns two polynomials say G and H so that G + H * alpha
  36 approximates the L-series where alpha is a root of x^2 - a_p x + p.
  37 
  38 The other main procedure is iwasawa_invariants_w_curve.  This
  39 takes 6 arguments, say N,V,p,a_p,n,D.
  40 
  41 Here N is the conductor of the curve
  42      V is the space of modular symbols of the curve
  43      p is some prime
  44      a_p is the a_p value of the curve
  45      n is the level of accuracy
  46      D represents twisting by a quadratic character of conductor D
  47 
  48 This procedure is supposed to return the iwasawa invariants of 
  49 the L-series.  
  50 
  51 Precisely, in the ordinary case it returns a 6-tuple of data,
  52 namely the rank of the curve, 
  53        the p-adic valuation of the leading non-zero term of the L-series,
  54        the mu invariant of the L-series
  55        the lambda invariant of the L-series
  56        the newton polygon of the L-series.
  57                  -this is returned as [ [m_1,s_1] , [m_2,s_2] , ... ]
  58                  where m_1 is the number of roots of slope s_1...;
  59                  a slope of 10000 represents the root 0.
  60        a boolean saying whether the above data is correct.
  61          (note if mu is positive then a "true" only means that the data
  62          is as accurate as it can be using a polynomial of degree
  63          p^(n-1)-1; that is it could have a another root of smaller slope)
  64 
  65 
  66 In the supersingular case, the p-adic L-series is not an iwasawa
  67 function so the above data couldn't make any sense.  Instead it
  68 returns an 11-tuple of data consisting of the rank and then the
  69 same data for the g and h functions i construct in my thesis.  i should
  70 have a rough draft of my thesis ready soon if you are interested
  71 in how they are defined.  basically, they are the "real" and "imaginary"
  72 parts of the L-series with their trivial zeroes removed.  They are
  73 in fact iwasawa functions (actually integral power series).
  74 
  75 I've been back and forth between NY and Boston for the last
  76 few weeks.  Now I should be back in Boston full time for the
  77 next month.
  78 
  79 Okay...let me try and send out the code.
  80 
  81 Robert
  82 *********************************************************************/
  83 
  84 R<x>:=PolynomialRing(Rationals());
  85 R<T>:=PolynomialRing(Rationals());
  86 
  87 //This function cuts out the subspace of H_1(X_0(N),Z) where
  88 //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
  89 //Input:  	N - conductor
  90 //		a_2 - eigenvalue of T_2
  91 //		a_3 - eigenvalue of T_3
  92 //		a_5 - eigenvalue of T_5
  93 //		a_7 - eigenvalue of T_7
  94 //
  95 //Output:	space of modular symbols of level N and specified eigenvalues
  96 function form_curve(N,a_2,a_3,a_5,a_7)
  97 
  98    R<x>:=PolynomialRing(Rationals());
  99    M:=ModularSymbols(N,2,+1);
 100    I:=[<2,x-a_2>,<3,x-a_3>,<5,x-a_5>,<7,x-a_7>]; 
 101    V:=Kernel(I,M);
 102    ModularSymbolEven(V,1/2);
 103 
 104    return(V);
 105 end function;
 106 
 107 function table_of_modular_symbols(V,p,n,D) 
 108    local q,table,temp,k,scale,q1,mid,temp2,e;
 109 
 110    if (D mod p eq 0) then
 111       print "Twisting by something not relatively prime to p!";
 112    end if;
 113 
 114    if (KroneckerSymbol(D,-1) eq -1) then
 115       print "Twisting by an odd character!";
 116    end if;
 117 
 118    q1 := p^(n-1);
 119    q := q1*p;
 120 
 121    k:=1;
 122    while (ModularSymbolEven(V,k/p^10)[1] eq 0) do
 123       k := k+1;
 124    end while;
 125 
 126    scale := Denominator(ModularSymbolEven(V,k/q)[1]);
 127    table := [];
 128 
 129    for a in [0..(q-1) div 2] do
 130       if (Valuation(a,p) le 2) or (a eq 0) then
 131          temp2 :=0;
 132          for b in [1..D] do
 133             e := KroneckerSymbol(D,b);
 134             if (e ne 0) then
 135                temp := ModularSymbolEven(V,(D*a-q*b)/(D*q))[1] * scale;
 136                if Denominator(temp) ne 1 then   
 137                   print "PROBLEMS!  There are still denominators in ",
 138                         "the modsyms.";
 139                end if;   
 140                temp2 := temp2 + KroneckerSymbol(D,b)*temp;
 141             end if;
 142          end for;
 143          table := table cat [Numerator(temp2)];
 144       else
 145          table := table cat [0];
 146       end if;
 147    end for;
 148 
 149    //USING the fact that [r/s]^+ equals [-r/s]^+
 150    mid := ((q-1) div 2) + 1;
 151    for a in [mid..q-1] do
 152       table := table cat [table[q-a+1]];
 153    end for;
 154    table := table cat [scale];
 155 
 156    return(table);
 157 end function;
 158 
 159 //Multiplies a and b interpreted as a[1] + sqrt(-p)*a[2].
 160 function multiply(a,b,p)
 161    return([a[1]*b[1] - p*a[2]*b[2],a[1]*b[2]+a[2]*b[1]]);
 162 end function;
 163 
 164 //Here answer[a] = [alpha^-a]
 165 //n+1 is the maxl power and m is the accuracy
 166 function powers_of_alpha(p,a_p,n,m)
 167    local temp,answer,t,q;
 168 
 169    if a_p mod p ne 0 then
 170       temp := a_p mod p;
 171 
 172       for k in [2..m] do
 173          t := -(temp^2-a_p*temp+p) div p^(k-1) mod p *
 174             InverseMod(2*temp - a_p,p);
 175          temp := temp + t*p^(k-1);
 176       end for;
 177 
 178       q := p^m;
 179       temp := InverseMod(temp,q);
 180       answer := [[temp]];
 181 
 182       for k in [2..n+1] do
 183          Append(~answer,[answer[k-1][1]*answer[1][1] mod q]);
 184       end for;
 185 
 186       return(answer);
 187    else
 188       if a_p eq 0 then
 189          temp := [0,-1/p];
 190       else
 191          if a_p eq -3 then
 192             temp := [-1/2,-1/6];
 193          else
 194             temp := [1/2,-1/6];
 195          end if;
 196       end if;
 197       answer := [temp];
 198       for k in [2..n+1] do
 199          Append(~answer,multiply(answer[k-1],answer[1],p));
 200       end for;
 201       return(answer);
 202    end if;
 203 end function;
 204 
 205 //Here answer = {a} mod p^m
 206 function one_teich(a,p,m)
 207    local answer,t;
 208 
 209    answer := a mod p;
 210    for k in [2..m] do
 211       t := 0;
 212       while (answer + t*p^(k-1))^(p-1) mod p^k ne 1 do
 213          t := t+1;
 214       end while;
 215       answer := answer + t*p^(k-1); 
 216    end for;
 217 
 218    return(answer);
 219 end function;
 220          
 221 //Here answer[a] = {a} (mod p^m) 
 222 function list_of_teich(p,m)
 223    local answer;
 224 
 225    answer := [];
 226    for a in [1..p-1] do
 227       Append(~answer,one_teich(a,p,m));
 228    end for;
 229    return(answer);
 230 end function;
 231 
 232 //Here answer[a] = gamma^(a-1)
 233 //SO THERE IS A SHIFT
 234 //n is the maxl power and m is the accuracy
 235 function powers_of_gamma(gamma,p,n)
 236    local answer;
 237 
 238    answer := [1];
 239    for k in [2..p^(n-1)] do
 240       Append(~answer,answer[k-1]*gamma mod (p^n));
 241    end for;
 242 
 243    return answer;
 244 end function;
 245 
 246 
 247 function lift_polynomial(f)
 248    return PolynomialRing(Rationals()) ! 
 249                 (PolynomialRing(Integers()) ! f);
 250 end function;
 251 
 252 //Here answer[a]=(1+T)^(a-1)
 253 //SO THERE IS A SHIFT
 254 //n is the maxl power and m is the accuracy
 255 function powers_of_oneplusT(p,n,m)
 256    local R,S,answer,q;
 257 
 258    q := p^m;
 259    R := ResidueClassRing(q);
 260    S<T> := PolynomialRing(R);
 261    answer := [1];
 262    for k in [2..p^(n-1)] do
 263       answer := answer cat [answer[k-1]*(1+T)] ;
 264    end for;
 265 
 266    return(answer);
 267 end function;
 268 
 269 function alpha(n,ALPHA)
 270      return(ALPHA[n]);
 271 end function;
 272 
 273 function modsym(a,MODSYM)
 274    return(MODSYM[a+1]);
 275 end function;
 276 
 277 function oneplusT(j,POWER)
 278    return(POWER[j+1]);
 279 end function;
 280 
 281 function teich(a,TEICH)
 282    return(TEICH[a]);
 283 end function;
 284 
 285 function gamma(j,GAMMA)
 286    return(GAMMA[j+1]);
 287 end function;
 288 
 289 function scale(list,c)
 290    local k,answer;
 291 
 292    answer:=[];
 293    for k in [1..#list] do
 294       answer := answer cat [list[k]*c];
 295    end for;
 296    return(answer);
 297 end function;
 298 
 299 function add(list1,list2)
 300    local k,answer;
 301    
 302    if #list1 ne #list2 then
 303       print "Error in add";
 304    end if;
 305 
 306    answer:=[];
 307    for k in [1..#list1] do
 308       answer := answer cat [list1[k]+list2[k]];
 309     end for;
 310 
 311    return(answer);
 312 end function;
 313 
 314 function subtract(list1,list2)
 315    local k,answer;
 316    
 317    if #list1 ne #list2 then
 318       print "Error in subtract";
 319    end if;
 320 
 321    answer:=[];
 322    for k in [1..#list1] do
 323       answer := answer cat [list1[k]-list2[k]];
 324     end for;
 325 
 326    return(answer);
 327 end function;
 328 
 329 function measure(a,p,n,ALPHA,MODSYM)
 330    local q,q1;
 331 
 332    q := p^n; 
 333    q1 := q div p;
 334    return( subtract(scale(alpha(n,ALPHA),modsym(a,MODSYM)),
 335            scale(alpha(n+1,ALPHA),modsym((a mod q1) * p ,MODSYM))));
 336 
 337 end function;
 338 
 339 //This procedure returns a polynomial approximation to the p-adic
 340 //L-series of the elliptic curve corresponding to V twisted by a quadratic
 341 //character of discriminant D.  The degree of the polynomial is p^n-1. 
 342 //If a_p is prime to p (ordinary) then the answer is simply a polynomial.
 343 //Otherwise (in the supersingular case) the answer is in the form [G,H]
 344 //where the L-series is approximated by G + sqrt(-p) * H.
 345 //
 346 //Input:   V - space of modular symbols
 347 //      p - prime
 348 //      a_p - eigenvalue of T_p acting on V
 349 //      n - relates to degree of approximation
 350 //      D - value to twist by
 351 function Lseries(V,p,a_p,n,D)
 352    local answer,q;
 353 
 354    ALPHA := powers_of_alpha(p,a_p,n,n+1);
 355    MODSYM := table_of_modular_symbols(V,p,n,D);
 356    POWER := powers_of_oneplusT(p,n,n+1);
 357    TEICH := list_of_teich(p,n+1);
 358    GAMMA := powers_of_gamma(1+p,p,n+1); 
 359 
 360    q1 := p^(n-1);
 361    q := q1*p;
 362    if a_p mod p ne 0 then
 363       answer := [0];
 364    else
 365       answer := [0,0];
 366    end if;
 367 
 368    for a in [1..p-1] do
 369       for j in [0..q1-1] do
 370          answer := add(answer, 
 371                    scale(measure(teich(a,TEICH)*gamma(j,GAMMA) mod 
 372                                 q,p,n,ALPHA,MODSYM),
 373                         lift_polynomial(oneplusT(j,POWER))));
 374       end for;
 375    end for;
 376 
 377    answer:=scale(answer,1/D);
 378    answer:=scale(answer,1/MODSYM[#MODSYM]);
 379 
 380    return(answer);
 381 end function;
 382 
 383 function highest_power(n,p)
 384    local k,temp;
 385 
 386    k:=0;
 387    temp:=1;
 388    while (n ge temp) do
 389       temp := temp*p;
 390       k := k+1;
 391    end while;
 392    return(k-1);
 393 end function;
 394 
 395 function valuation_of_coefficients(f,p,n)
 396    local answer,j,temp,okay,bound,list_bad;
 397 
 398    answer := [];
 399    okay:=true;
 400    list_bad:=[];
 401    for j in [1..Degree(f)+1] do
 402       if Coefficient(f,j-1) ne 0 then
 403          answer := answer cat [Valuation(Coefficient(f,j-1),p)];
 404       else 
 405          answer := answer cat [10000];
 406       end if;
 407       
 408       bound := n-1-highest_power(j-1,p);
 409       if (Valuation(Coefficient(f,j-1),p) ge bound) and (j ne 1) then
 410          list_bad := list_bad cat [j];
 411       end if;
 412    end for;
 413 
 414    if f eq 0 then
 415       answer := [];
 416       answer := answer cat [10000] cat [10000] cat [10000];
 417    end if;
 418    
 419    return(<answer,list_bad>);
 420 end function;
 421    
 422 function newton_polygon(C)
 423    local k,here,next,min_slope,answer,okay,prev_slope;
 424 
 425    answer:= [];
 426    okay:=true;
 427    here:=1;
 428    next:=1;
 429 
 430    while here lt #C do
 431       min_slope:=100000;
 432       for k in [here+1..#C] do
 433          if (C[k]-C[here])*(1.0)/(k-here) le 1.0*min_slope then
 434             min_slope := (C[k]-C[here])/(k-here);
 435             next := k;
 436          end if;
 437       end for;
 438       answer := answer cat [[next-here,-min_slope]];
 439       here := next;
 440    end while;   
 441    return(answer);
 442 end function;
 443 
 444 function find_last_nontriv(NP,p,sign)
 445    local slope,place;
 446 
 447    place := 0;
 448    if sign mod 2 eq 0 then
 449       slope := 1/(p*(p-1));
 450    else
 451       slope := 1/(p-1);
 452    end if;
 453 
 454    for k in [1..#NP] do
 455       if (NP[k][2] gt 0) then
 456          if (NP[k][2] ne slope) or (NP[k][1] ne 1/slope) then
 457             place:=k;
 458          else
 459             slope:=slope/(p^2);
 460          end if;
 461       end if;
 462    end for;
 463 
 464    return(place);
 465 end function;
 466 
 467 
 468 function refine_newton_polygon(NP,p,a_p,sign)
 469    local answer,k,last;
 470    
 471    answer := [];
 472    if a_p mod p ne 0 then
 473       for k in [1..#NP] do
 474          if NP[k][2] gt 0 then
 475             answer := answer cat [NP[k]];
 476          end if;
 477       end for;
 478    else
 479       last := find_last_nontriv(NP,p,sign);
 480       for k in [1..last] do
 481          answer := answer cat [NP[k]];
 482       end for;
 483    end if;
 484             
 485    return(answer);
 486 end function;
 487 
 488 function minimize_ord(C,p,n) 
 489    local k,start;
 490 
 491    start:=1;
 492    if C[1] eq 10000 then
 493       start := 2;
 494    end if;
 495 
 496    for k in [start..#C] do
 497       if C[k] ge n-1-highest_power(k-1,p) then
 498          C[k] := n-1-highest_power(k-1,p);
 499       end if;
 500    end for;
 501    return(C);
 502 end function;
 503 
 504 function in_list(list,num)
 505    local k,flag;
 506 
 507    flag:=false;
 508    for k in [1..#list] do
 509       if list[k] eq num then
 510          flag := true;
 511       end if;
 512    end for;
 513    
 514    return(flag);
 515 end function;
 516 
 517 function compute_rank(C,N,r,D)
 518    local rank;
 519 
 520    if C[1] eq 10000 then
 521       rank := 1;
 522    else
 523       rank :=0;
 524     end if;
 525 
 526    if (r eq 0) and (rank eq 1) and (KroneckerSymbol(D,N) eq 1) then  
 527           //it has even rank >= 2.
 528       rank:=2;
 529       C[2]:=10000;
 530    end if;
 531 
 532    if (r eq 1) and (rank eq 1) and (KroneckerSymbol(D,N) eq -1) then  
 533           //it has even rank >= 2.
 534       rank:=2;
 535       C[2]:=10000;
 536    end if;
 537 
 538    return(<C,rank>);
 539 end function;
 540 
 541 
 542 function check_NP(NP,C,list_bad,rank,p,a_p,n,sign)
 543    local temp,val,k,good;
 544 
 545    //checks if break points are accurate.
 546    
 547    good:=true;
 548    temp := refine_newton_polygon(NP,p,a_p,sign);
 549    val := 1;
 550    for k in [1..#temp] do
 551       val := val + temp[k][1];
 552       if (in_list(list_bad,val)) then
 553          good := false;
 554          print "Failed at a breakpoint,",val;
 555       end if;
 556    end for;
 557 
 558    //checks if lowering the valuation of any unknown coefficient to its
 559    //minimal possible value affects the newton polygon.
 560 
 561    C_min := minimize_ord(C,p,n); //lowers all unknown coefficients 
 562                                  // to lowest possible valuation
 563    
 564    if rank eq 2 then
 565          C_min[2] := 10000;
 566    end if;
 567 
 568    NP_min := newton_polygon(C_min);
 569    if refine_newton_polygon(NP_min,p,a_p,sign) 
 570           ne refine_newton_polygon(NP,p,a_p,sign) then
 571       good := false;
 572       print "Failed while lowering.";   
 573    end if;
 574 
 575    return(good);
 576 end function;
 577 
 578 function lambda(NP)
 579    local k,lam;
 580 
 581    k:=1;
 582    lam := 0;
 583    while (k le #NP) and (NP[k][2] gt 0) do
 584         lam := lam + NP[k][1];
 585       k := k+1;
 586    end while;
 587 
 588    return(lam);
 589 end function;
 590 
 591 function mu(NP,C,rank)
 592    local k,val,m,start;
 593 
 594    k:=1;
 595    val := 0;
 596    if (rank ne 2) then
 597       while (k le #NP) and (NP[k][2] gt 0) do
 598          if NP[k][2] lt 1000 then
 599             val := val + NP[k][1]*NP[k][2];
 600          end if;
 601          k := k+1;
 602       end while;
 603    else
 604       if (#NP gt 0) and (NP[1][1] eq 1) then 
 605          start := 3;
 606       else
 607          if #NP eq 0 then
 608             start := 1;
 609          else
 610             start := 2;
 611          end if;
 612       end if;
 613       k:=start;
 614       while (k le #NP) and (NP[k][2] gt 0) do
 615          if NP[k][2] lt 1000 then
 616             val := val + NP[k][1]*NP[k][2];
 617          end if;
 618          k := k+1;
 619       end while;
 620    end if;
 621 
 622    if rank eq 0 then
 623       m := C[1]-val;
 624    else
 625       if rank eq 1 then
 626          m := C[2] - val;
 627       else
 628          m := C[3] - val;
 629       end if;
 630    end if;
 631 
 632    return(m);
 633 
 634 end function;
 635 
 636 function remove_two(list)
 637    local answer,k;
 638 
 639    answer := [];
 640 
 641    for k in [1..#list] do
 642       if list[k] ne 2 then
 643          answer := answer cat [list[k]];
 644       end if;
 645    end for;
 646 
 647    return(answer);
 648 
 649 end function;
 650 
 651 function leading_term(C,rank)
 652 
 653    if rank eq 0 then
 654       return(C[1]);
 655    else
 656       if rank eq 1 then
 657          return(C[2]);
 658       else
 659          if rank eq 2 then
 660             return(C[3]);
 661          end if;
 662       end if;
 663    end if;
 664 end function;
 665 
 666 function compute_untwisted_rank(V,p,a_p)
 667    local L,C;
 668 
 669    L:=Lseries(V,p,a_p,2,1);
 670    L:=L[1];
 671    C := valuation_of_coefficients(L,p,1);
 672    if C[1][1] ne 10000 then
 673       return(0);
 674    else 
 675       return(1);
 676    end if;
 677 end function;
 678 
 679 function remove_unit_zeroes(NP)
 680    local k,answer;
 681 
 682    answer := [];
 683    for k in [1..#NP] do
 684       if NP[k][2] gt 0 then
 685          answer := answer cat [NP[k]];
 686       end if;
 687    end for;
 688    return(answer);
 689 end function;
 690 
 691 //Has the space of modular symbols already pre-computed
 692 function iwasawa_invariants_w_curve_ord(N,V,p,a_p,n,D)
 693    local C,NP,k,lam,m,rank,j,list_bad,list_breapoints,
 694          C_min,NP_min,val,good,m2,lam2,temp,leading;
 695 
 696    L:=Lseries(V,p,a_p,n,D)[1];
 697    C := valuation_of_coefficients(L,p,n);
 698    list_bad:=C[2];
 699    C:=C[1];
 700       
 701    //FIRST COMPUTE RANK
 702    
 703    rank := compute_untwisted_rank(V,p,a_p);
 704    C := compute_rank(C,N,rank,D);
 705 
 706    rank := C[2];  //changes valuation of a_1 to 10000 if guesses rank 2.
 707    C := C[1];   
 708 
 709    if rank eq 2 then
 710       list_bad := remove_two(list_bad);
 711    end if;
 712 
 713    NP := newton_polygon(C);
 714    good := check_NP(NP,C,list_bad,rank,p,a_p,n,0);  //0 represents nothing.
 715    NP := remove_unit_zeroes(NP);
 716    leading := leading_term(C,rank);
 717    
 718    if (leading eq 0) and (in_list(list_bad,rank+1) eq false) then
 719         good:=true;
 720    end if;
 721 
 722    lam:=lambda(NP);
 723    m := mu(NP,C,rank);   
 724 
 725    return(<rank,leading,m,lam,NP,good>);
 726 end function;
 727 
 728 function scale_by_p(C)
 729    local k;
 730 
 731    for k in [1..#C] do
 732       if C[k] ne 10000 then
 733          C[k]:=C[k]+1;
 734       end if;
 735    end for;
 736 
 737    return(C);
 738 end function;
 739 
 740    
 741 function remove_forced_zeroes(NP,p,init_slope)
 742    local k,slope,temp;
 743 
 744    temp := [];
 745    k := 1;   
 746    slope:=init_slope;
 747    while (k le #NP) and (NP[k][2] gt 0)  do
 748       if NP[k][2] ne slope then
 749          temp := temp cat [NP[k]];
 750       else
 751          if NP[k][1] - 1/slope ne 0 then
 752             temp := temp cat [[NP[k][1] - 1/slope,slope]];
 753          end if;
 754           slope := slope/(p^2);
 755       end if;
 756       k := k+1;
 757    end while;
 758 
 759    return(temp);
 760 end function;
 761 
 762 //Has the space of modular symbols already pre-computed
 763 function iwasawa_invariants_w_curve_ss(N,V,p,a_p,n,D)
 764    local L,G,H,k,CG,CH,mg,mh,lamg,lamh,NPG,NPH,rank,rankg,rankh,
 765       list_badg,list_badh,denomg,denomh,goodg,goodh,
 766       list_breakpointsg,list_breakpointsh,
 767       leadingg,leadingh;
 768 
 769    L := Lseries(V,p,a_p,n,D);
 770    G := L[1];
 771    H := L[2];
 772 
 773    //   denomg and denomh represent the number of p's 
 774    //   in the denominator of the measure
 775 
 776    if n mod 2 eq 0 then
 777       denomg:=n div 2;
 778       denomh:=(n+2) div 2;
 779    else
 780       denomg:=(n+1) div 2;
 781       denomh:=(n+1) div 2;
 782    end if;   
 783       
 784    CG := valuation_of_coefficients(G,p,n-denomg);
 785    CH := valuation_of_coefficients(H,p,n-denomh);
 786 
 787    list_badg := CG[2];
 788    list_badh := CH[2];
 789 
 790    CG:=CG[1];
 791    CH:=CH[1];
 792 
 793    //Computing the rank
 794 
 795    rank:=compute_untwisted_rank(V,p,a_p);
 796 
 797    CG := compute_rank(CG,N,rank,D); //changes a_1 if guesses rank 2
 798    rankg := CG[2];
 799    CG := CG[1];
 800 
 801    CH := compute_rank(CH,N,rank,D);
 802    rankh := CH[2];
 803    CH := CH[1];
 804 
 805    if rankg ne rankh then
 806       print "WEIRD..ranks different in G and H";
 807       PrintFile("notes","ranks different in G and H");
 808       PrintFile("notes",D);
 809    end if;
 810 
 811    rank := rankg;
 812    if rank eq 2 then
 813       list_badg := remove_two(list_badg);
 814    end if;
 815 
 816    if rank eq 2 then
 817       list_badh := remove_two(list_badh);
 818    end if;
 819 
 820    NPG := newton_polygon(CG);
 821    NPH := newton_polygon(CH);
 822 
 823    goodg := check_NP(NPG,CG,list_badg,rank,p,a_p,n,0); //0 represents g
 824    goodh := check_NP(NPH,CH,list_badh,rank,p,a_p,n,1); //1 represents h
 825 
 826    leadingg := leading_term(CG,rank);
 827    leadingh := leading_term(CH,rank);
 828 
 829    if (a_p eq 0) then
 830       if (leadingg eq -1) and (in_list(list_badg,rank+1) eq false) then
 831           goodg:=true;
 832       end if;
 833 
 834       if (leadingh eq -1) and (in_list(list_badh,rank+1) eq false) then
 835           goodh:=true;
 836       end if;
 837    end if;
 838 
 839 
 840    if a_p eq 0 then
 841       NPG := remove_forced_zeroes(NPG,p,1/((p-1)*p));
 842       NPH := remove_forced_zeroes(NPH,p,1/(p-1));
 843    else
 844       if a_p eq -3 then
 845          NPG := remove_forced_zeroes(NPG,p,1/((p-1)*p));
 846          NPH := remove_forced_zeroes(NPH,p,1/(p-1)*p^2);
 847       end if;
 848    end if;
 849 
 850    lamg := lambda(NPG);
 851    lamh := lambda(NPH);   
 852 
 853    mg := mu(NPG,CG,rank)+1;
 854    mh := mu(NPH,CH,rank)+1;
 855 
 856    return(<rank,leadingg,mg,lamg,NPG,goodg,leadingh,mh,lamh,NPH,goodh>);
 857 end function;
 858 
 859 // This procedure returns many invariants of the curve corresponding
 860 // to V twisted by D.  
 861 //
 862 // In the ordinary case the data returned is:
 863 //   <rank,valuation of leading coefficient,mu invariant,lambda invariant,
 864 //    newton polygon of L-series, boolean saying whether the above data
 865 //    are accurate>
 866 //
 867 // The newton polygon is returned as [n_1,s_1],[n_2,s_2],... where there are
 868 // n_1 roots of slope s_1, etc.
 869 //
 870 // In the supersingular case the data returned is similar.  The first entry is
 871 // the rank, then the next five are the data for my "g" iwasawa function and
 872 // then the next five for my "h" iwasawa function.
 873 
 874 function iwasawa_invariants_w_curve(N,V,p,a_p,n,D)
 875    if a_p mod p ne 0 then
 876       return(iwasawa_invariants_w_curve_ord(N,V,p,a_p,n,D));
 877    else
 878       return(iwasawa_invariants_w_curve_ss(N,V,p,a_p,n,D));
 879    end if;
 880 end function;