CoCalc Public Fileswww / cgi-bin / mfd / buhler_gross.gp
Author: William A. Stein
Compute Environment: Ubuntu 18.04 (Deprecated)
1\\ Buhler-Gross algorithm for computing L-series derivatives
2
3\\ Implementation by [email protected]
4\\ Repaired and improved by [email protected]
5
6\\ This is practical for any rank and conductors <= 10^10 in 10
7\\ minutes per L' calculation (on P3/500E).
8
9default(primelimit,2000000);
10
11\\ Section 0: prime-number utility which really should be
12\\ built-in. This version due to Karim (reference for upper and lower
13\\ bounds is Rosser & Schoenfeld, Approximate formulas for some functions
14\\ of prime numbers, Illinois J Math 6 (1962) pp 64-94)
15
16pith(x) =
17{
18  local(l,r,lx);
19  if (x <= 2, return (x==2));
20  lx = log(x);
21  l = floor(x/lx);
22  r = floor(l * (1 + 3/(2*lx)));
23  while (r-l>1,
24    m = (l+r)>>1;
25    if (prime(m)<=x, l=m, r=m));
26  l;
27}
28
29\\ Section 1: G_r
30
31\\ We need a load of Magic Constants, which are the coefficients
32\\ of the Taylor expansion of GAMMA(1+s). Now, [BGZ] we have
33\\ log(Gamma(1+s)) = -gamma s + sum(n=2,infinity,(-1)^n*zeta(n)/n)
34\\ and Pari can do polynomial exponential
35
36default(seriesprecision,30);
37ellld_ll = -Euler*s + sum(n=2,30,(-1)^n * zeta(n) / n * s^n);
38ellld_mm = exp(ellld_ll);
39ellld_pn_coeffs=vector(25,i,polcoeff(ellld_mm,i-1));
40
41ellld_P(r,t)=local(m);sum(m=0,r,ellld_pn_coeffs[r-m+1]*t^m/factorial(m));
42
43{ellld_G(r,x)=
44   local(ss,oldss,n,sn);
45   if (x>30,return(0));
46   ss = ellld_P(r,log(1/x));
47   oldss = ss-1;n=1;sn=-((-1)^r);
48   while (abs(ss-oldss)>1e-15,
49     oldss = ss;
50     ss = ss + sn*x^n/(n^r*factorial(n));
51     n=n+1;sn=-sn;
52  );
53  2*ss;
54}
55
56\\ Section 2: the BG recursion
57\\ Some global variables are necessary, but I've named them
58\\ to try to avoid contaminating the namespace *too* much.
59
60ellld_getap(r)=if(r<=ellld_pnum,ellld_ap[r],ellap(est,prime(r)))
61ellld_getp(r)=if(r<=ellld_pnum,ellld_p[r],prime(r))
62
63{
65	local(j,j0,next_a);
66	if(a==0,j0=i,ellld_sum=ellld_sum+a*ellld_G(ellld_r,ellld_X*n)/n;j0=1;);
67	j=j0;
68	while(j<=i && ellld_getp(j)*n <= ellld_bnd,
69		next_a=a*ellld_getap(j);
70		if(j==i && ellld_N%ellld_getp(j)!=0,
71			next_a = next_a - ellld_getp(j)*last_a);
73		j=j+1;
74	)
75}
76
77{
78elllderiv(e,r,bnd) =
79  local(cached_pmax,loop_p_num,ithresh,N,est);
80  est	       = ellinit(e);
81  cached_pmax  = sqrt(bnd);
82  ellld_pnum   = pith(cached_pmax);
83  loop_p_num   = pith(bnd);
84  ellld_p      = primes(ellld_pnum);
85  ellld_ap     = vector(ellld_pnum,i,ellap(est,ellld_p[i]));
86\\ Version 2.0.19 has acknowledged bug in ellap for p=2
87  ellld_ap[1]  = ellak(est,2);
88  N	       = ellglobalred(e)[1];
89
90  ellld_N      = N;
91  ellld_X      = 2*Pi/sqrt(N);
92  ellld_r      = r;
93  ellld_bnd    = bnd;
94  ellld_sum    = ellld_G(r,ellld_X);
95  ithresh=2;
96  for(i=1,loop_p_num,
97\\    if(i==ithresh,
98\\      print([i,ellld_sum]);
99\\      ithresh=ithresh*2;
100\\    );
102  );
103  ellld_sum;
104}
105
106{
107  ellanalyticrank(e)=
108  local(rksgn,ZB,egrd,cond,localbnd,Lr1,RS,om);
109  egrd=ellglobalred(e);
110\\ We had to compute this anyway to get the conductor, and
111\\ ellap gives wrong answers on non-minimal curves ...
112  if(egrd[2]!=[1,0,0,0],
113    print("INPUT CURVE NOT MINIMAL");
114    e=ellchangecurve(e,egrd[2]);
115    print("MINIMISED TO ",[e[1],e[2],e[3],e[4],e[5]]));
116  cond=egrd[1];
117\\ No, I can't justify this, but it should be big enough ...
118  localbnd=round(2*sqrt(cond));
119  print("Summing ",localbnd," a_n terms");
120  rksgn=ellrootno(e);
121  if(rksgn==-1,print("Rank is odd"),print("Rank is even"));
122  ZB=(1-rksgn)/2;
123  Lr1=0;
124\\ And I'm sure it's conjectural that this is close-enough to 0
125  while(abs(Lr1)<0.001,
126    Lr1=elllderiv(e,ZB,localbnd);
127    print("L^(",ZB,")=",Lr1);
128    ZB=ZB+2);
129\\ Also compute the quantity reg*sha from the BSD conjecture
130  om=e[15]; if(e.disc>0,om=2*om);
131  [ZB-2,Lr1,(Lr1*elltors(e)[1]^2)/(om*egrd[3])];
132}
133
134{ellRegSha(e,r)=
135  local(egrd,om,Lfn);
136  om=e[15];if(e.disc>0,om=2*om);
137  egrd=ellglobalred(e);
138  Lfn=elllderiv(e,r,2*sqrt(egrd[1]));
139  (Lfn*elltors(e)[1]^2)/(om*egrd[3]);
140}
141
142ellSha(e,pts)=ellRegSha(e,length(pts))/matdet(ellheightmatrix(e,pts));
143
144
145
146