CoCalc Shared Fileswww / cgi-bin / mfd / buhler_gross.gpOpen in CoCalc with one click!
Author: William A. Stein
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
9
default(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
16
pith(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
36
default(seriesprecision,30);
37
ellld_ll = -Euler*s + sum(n=2,30,(-1)^n * zeta(n) / n * s^n);
38
ellld_mm = exp(ellld_ll);
39
ellld_pn_coeffs=vector(25,i,polcoeff(ellld_mm,i-1));
40
41
ellld_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
60
ellld_getap(r)=if(r<=ellld_pnum,ellld_ap[r],ellap(est,prime(r)))
61
ellld_getp(r)=if(r<=ellld_pnum,ellld_p[r],prime(r))
62
63
{
64
BGadd(n,i,a,last_a) =
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);
72
BGadd(ellld_getp(j)*n,j,next_a,a);
73
j=j+1;
74
)
75
}
76
77
{
78
elllderiv(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
\\ );
101
BGadd(ellld_getp(i),i,ellld_getap(i),1)
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
142
ellSha(e,pts)=ellRegSha(e,length(pts))/matdet(ellheightmatrix(e,pts));
143
144
145
146