Open in CoCalc
1
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
2
\\
3
\\ The program modsym.gp
4
\\
5
\\
6
\\ HISTORY:
7
\\
8
\\ *) Original program written for PARI 1.x by Joseph L. Wetherell
9
\\
10
\\ *) July 1999: William Stein updated the program for PARI 2.x.
11
\\
12
\\ *) March 2001: Dominique Bernardi and Bernadette Perrin-Riou made
13
\\ the modifications described below:
14
\\
15
\\ Les modifications sont de type : dans le programme modsym.gp,
16
\\ les relations binaires n'�taient pas toutes prises en compte,
17
\\ par exemple, certaines relations du type 2*x = 0 (voir par
18
\\ exemple le cas N=6, N=106). Nous avons donc d� enlever
19
\\ l'astuce qui permettait de r�duire les matrices de taille
20
\\ N � une taille N/2, ce qui ralentit bien s�r le programme.
21
\\ Nous avons d'autre part rajout� un signe +1 ou -1 afin
22
\\ de calculer aussi l'application de l'espace des formes
23
\\ modulaires dans la partie + ou - de l'homologie.
24
\\ Enfin, nous avons rajout� deux proc�dures ellsym et
25
\\ xpm qui permettent de calculer les composantes x^\pm E
26
\\ du symbole modulaire associ�e � une courbe elliptique
27
\\ donn�e dans GP/Pari par ellinit ou plut�t � la forme
28
\\ modulaire asssoci�e. Le but est de calculer des valeurs
29
\\ de fonctions L p-adiques associ�es � une courbe elliptique
30
\\ et � un discriminant de corps quadratiques (positif ou n�gatif).
31
\\
32
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
33
34
talk=1; \\ verbosity flag.
35
36
{mnormal(cond,v)=
37
lift(v*Mod(1,cond));
38
}
39
40
{generatemsymbols(cond,
41
num,ret,curn)=
42
\\ num = [\Gamma:\Gamma_0(N)] = N * Prod_{p|cond} (1+p^-1)
43
num = cond;
44
v = factor(cond);
45
for(n=1,length(v~),
46
num = num*(1+1/v[n,1]);
47
);
48
49
\\ initialize M-symbol list
50
ret = vector(num,c,0);
51
curn = 1;
52
53
\\ generate M-symbols in three lists:
54
\\ list 1: (c:1) for 0 <= c < N
55
for(c=0, cond-1,
56
ret[curn] = [c,1];
57
curn = curn + 1;
58
);
59
60
\\ list 2: (1:d) for 0 <= d < N and gcd(d,N)>1
61
for(d=0, cond-1,
62
if(gcd(d,cond)>1,
63
ret[curn] = [1,d];
64
curn = curn + 1;
65
,);
66
);
67
68
\\ list 3: (c:d) with c|N, c<>1,N, gcd(c,d)=1, gcd(d,N)>1
69
v = curn;
70
for(d=2, cond-1,
71
if(gcd(d,cond)>1,
72
fordiv(cond,c,
73
if(c<>1 && c<>cond && gcd(c,d)==1,
74
if(0 == sum(n=v,curn-1, ((ret[n][1]*d-c*ret[n][2])%cond)==0 ),
75
ret[curn] = [c,d];
76
curn = curn + 1;
77
,);
78
,);
79
);
80
,);
81
);
82
83
if (curn -num - 1, print("pas le bon nombre de symboles"));
84
ret;
85
}
86
87
{hashmsymbols(cond,msymbols,
88
ret,tmp) =
89
ret = matrix(cond,cond,c1,d1,0);
90
for(n=1,cond-1,
91
if(gcd(n,cond)==1,
92
for(k=1,length(msymbols),
93
tmp = [msymbols[k][1]*n,msymbols[k][2]*n];
94
tmp = mnormal(cond,tmp);
95
ret[tmp[1]+1,tmp[2]+1] = k;
96
);
97
,);
98
);
99
ret;
100
}
101
102
{dohash(cond,hash,v,
103
tmp)=
104
v=mnormal(cond,v);
105
hash[v[1]+1,v[2]+1];
106
}
107
108
doV(v) = [-v[1],v[2]]
109
doS(v) = [-v[2],v[1]]
110
doT1(v) = [-v[1]-v[2],v[1]]
111
doT2(v) = [v[2],-v[2]-v[1]]
112
113
{msymbolrelations(cond,msymbols,hash,signe,
114
mask1,n1,n2,n3,s1,s2,s3,relat,nsym,freebasis,symtofree,
115
tempr) =
116
nsym = length(msymbols);
117
relat = matid(nsym);
118
mask1 = vector(nsym,n,1);
119
120
121
for(n=1,nsym,
122
n1 = msymbols[n];
123
n2 = doS(msymbols[n]);
124
n1 = dohash(cond,hash,n1);
125
n2 = dohash(cond,hash,n2);
126
r = relat[n1,] + relat[n2,];
127
minc = vecmax(abs(r));
128
for(k=1,nsym,
129
if(r[k] && abs(minc)>=abs(r[k]),
130
minc = r[k];
131
mincloc = k;
132
,);
133
);
134
if(minc,
135
mask1[mincloc] = 0;
136
r = r/minc;
137
relat = relat - relat[,mincloc]*r;
138
,);
139
);
140
141
for(n=1,nsym,
142
n1 = msymbols[n];
143
n2 = doT1(msymbols[n]);
144
n3 = doT2(msymbols[n]);
145
n1 = dohash(cond,hash,n1);
146
n2 = dohash(cond,hash,n2);
147
n3 = dohash(cond,hash,n3);
148
r = relat[n1,] + relat[n2,] + relat[n3,];
149
minc = vecmax(abs(r));
150
for(k=1,nsym,
151
if(r[k] && abs(minc)>=abs(r[k]),
152
minc = r[k];
153
mincloc = k;
154
,);
155
);
156
if(minc,
157
mask1[mincloc] = 0;
158
r = r/minc;
159
relat = relat - relat[,mincloc]*r;
160
,);
161
);
162
163
for(n=1,nsym,
164
n1 = msymbols[n];
165
n2 = doV(msymbols[n]);
166
n1 = dohash(cond,hash,n1);
167
n2 = dohash(cond,hash,n2);
168
r = relat[n1,] -signe*relat[n2,];
169
minc = vecmax(abs(r));
170
for(k=1,nsym,
171
if(r[k] && abs(minc)>=abs(r[k]),
172
minc = r[k];
173
mincloc = k;
174
,);
175
);
176
if(minc,
177
mask1[mincloc] = 0;
178
r = r/minc;
179
relat = relat - relat[,mincloc]*r;
180
,);
181
);
182
183
n1 = vector(nsym,n,mask1[n]*n);
184
n2 = eval(Set(simplify(n1)));
185
freebasis = vector(length(n2)-1,n,n2[n+1]);
186
symtofree = matrix(nsym,length(freebasis),n1,n2,relat[n1,freebasis[n2]]);
187
[freebasis,symtofree~];
188
}
189
190
cusps;
191
{initcusps() = cusps=[];}
192
193
{testcuspeq(cond,cusp1,cusp2,
194
p1,q1,p2,q2,s1,s2,n)=
195
p1 = cusp1[1];
196
q1 = cusp1[2];
197
p2 = cusp2[1];
198
q2 = cusp2[2];
199
s1 = if(q1>2, lift(Mod(1,q1)/Mod(p1,q1)), 1);
200
s2 = if(q2>2, lift(Mod(1,q2)/Mod(p2,q2)), 1);
201
202
foo=Mod(s1*q2-s2*q1, gcd(q1*q2,cond))==0;
203
204
foo;
205
}
206
207
{getcusp(cond,v,signe,
208
ret,negv) =
209
ret = 0;
210
negv = [-v[1],v[2]];
211
for(n=1,length(cusps),
212
if(testcuspeq(cond,v,cusps[n][1]), ret=[n,cusps[n][2]],
213
if(testcuspeq(cond,negv,cusps[n][1]), ret=[n,signe*cusps[n][2]]));
214
);
215
if(ret==0,
216
cusps = concat(cusps, [[v, 1-((signe < 0) && testcuspeq(cond,v,negv))]]);
217
ret = length(cusps);
218
ret = [ret,cusps[ret][2]];
219
,);
220
ret;
221
}
222
223
{delta(cond,msymbol,
224
v,c1,c2,ret)=
225
v = bezout(msymbol[1],msymbol[2]);
226
if(v[3]<>1, print("msymbol not coprime:", msymbol);1/0 ,);
227
c1 = [v[2],msymbol[1]];
228
c2 = [-v[1],msymbol[2]];
229
ret=[c2/gcd(c2[1],c2[2]), c1/gcd(c1[1],c1[2])];
230
ret;
231
}
232
233
{kerdelta(cond,msymbols,hash,freebasis,signe,
234
rels,trels,m,k,kinv,ksz)=
235
initcusps();
236
rels = vector(length(freebasis),n, delta(cond,msymbols[freebasis[n]]) );
237
for(n=1,length(rels),
238
rels[n] = [getcusp(cond,rels[n][1],signe),getcusp(cond,rels[n][2],signe)];
239
);
240
m = matrix(length(cusps),length(freebasis),n1,n2,0);
241
for(n=1,length(rels),
242
m[rels[n][1][1],n] = m[rels[n][1][1],n] + rels[n][1][2];
243
m[rels[n][2][1],n] = m[rels[n][2][1],n] - rels[n][2][2];
244
);
245
k = matker(m);
246
ksz = matsize(k);
247
kinvr = (matsupplement(k))^(-1);
248
[k, vecextract(kinvr,2^ksz[2]-1,2^ksz[1]-1)];
249
}
250
251
{modulartofree(cond,hash,symtofree,v,
252
q,hashval,n1,s1,sym,mods,ret)=
253
if(v[2]==0, v=[1,cond]);
254
q = contfrac(v[1]/v[2]);
255
q[1]=1;
256
for(n=3,length(q),
257
q[n]=q[n]*q[n-1]+q[n-2];
258
);
259
sym = vector(matsize(symtofree)[2],n,0)~;
260
for(n=2,length(q),
261
mods = [(-1)^(n)*q[n],q[n-1],1];
262
n1 = dohash(cond,hash,mods);
263
sym[n1] = sym[n1]+1;
264
);
265
ret=symtofree*sym;
266
ret;
267
}
268
269
{tpmatrix(cond,msymbols,hash,relations,hdata,p,
270
freebasis,symtofree,hbase,hdual,
271
sz,ret,col,sym,cusp,a)=
272
if(!isprime(p) || cond%p == 0,
273
print("Error, this program can only compute T_p for p prime and
274
coprime to N.");
275
return([]);
276
);
277
freebasis = relations[1];
278
symtofree = relations[2];
279
hbase = hdata[1];
280
hdual = hdata[2];
281
282
sz = matsize(hbase);
283
ret = matrix(sz[2],sz[2],n1,n2,0);
284
285
for(c=1,sz[2],
286
col = vector(sz[1],n,0)~;
287
for(r=1,sz[1],
288
289
a = hbase[r,c];
290
if(a,
291
sym = msymbols[freebasis[r]];
292
cusp = delta(cond,sym);
293
col = col -
294
a*modulartofree(cond,hash,symtofree,cusp[1]*[p,0;0,1]);
295
col = col +
296
a*modulartofree(cond,hash,symtofree,cusp[2]*[p,0;0,1]);
297
for(n=0,p-1,
298
col = col -
299
a*modulartofree(cond,hash,symtofree,cusp[1]*[1,0;n,p]);
300
col = col +
301
a*modulartofree(cond,hash,symtofree,cusp[2]*[1,0;n,p]);
302
);
303
,);
304
);
305
ret[,c] = hdual*col;
306
);
307
308
ret;
309
}
310
311
{gammap(cond,msymbols,hash,relations,hdata,p,
312
freebasis,symtofree,hbase,hdual,
313
sz,ret,col,sym,cusp,a)=
314
freebasis = relations[1];
315
symtofree = relations[2];
316
hbase = hdata[1];
317
hdual = hdata[2];
318
319
sz = matsize(hbase);
320
321
col = vector(sz[1],n,0)~;
322
for(n=1,p-1,
323
col = col + modulartofree(cond,hash,symtofree,[n,p]);
324
);
325
326
hdual*col;
327
}
328
329
defaultmodsymspace=0;
330
{modsym(N,signe,
331
symbols, hash, relations, hdata,
332
t1,t2,t3,t4)=
333
gettime();
334
if(talk,print1("1. Generating M-symbols \t\t("));
335
symbols = generatemsymbols(N);
336
if(talk,print(t1=gettime()," ms)");
337
print1("2. Hashing M-symbols \t\t\t("));
338
gettime();
339
hash = hashmsymbols(N,symbols);
340
if(talk,print(t2=gettime()," ms)");
341
print1("3. Quotienting out by relations \t("));
342
gettime();
343
relations = msymbolrelations(N,symbols,hash,signe);
344
if(talk,print(t3=gettime()," ms)");
345
print1("4. Computing the kernel of delta \t("));
346
gettime();
347
hdata = kerdelta(N,symbols,hash,relations[1],signe);
348
if(talk,print(t4=gettime()," ms)"));
349
if(talk,print("Total time.......................\t(",
350
t1+t2+t3+t4," ms)"));
351
352
defaultmodsymspace=[N,symbols,hash,relations,hdata,signe];
353
}
354
355
{T(p, M=0)=
356
if(type(M)=="t_INT",M=defaultmodsymspace;
357
if(M==0,print("T: Error -- no default space.")));
358
if(matsize(M)[2] != 6,
359
print("T: invalid modular symbols object."));
360
tpmatrix(M[1],M[2],M[3],M[4],M[5],p,M[6]);
361
}
362
363
{ellsym (ell,signe,
364
modele, cond, M, r, temp, res)=
365
modele = ellglobalred(ell);
366
cond = modele [1];
367
M = modsym (cond, signe);
368
r = matsize (M[5][1]) [2]; print("r = ", r);
369
res = if (r > 1, temp = [;];
370
forprime (p = 2, 10000,
371
if (matrank (temp) == r-1, break,);
372
print ("p = ", p);
373
if (cond % p, temp = concat (temp, T (p, M) - ellap(ell, p)*matid(r)),)
374
);
375
mattranspose (matker (mattranspose (temp))) * M[5][2],
376
M[5][2]
377
);
378
defaultellsyms = [cond, M[3], M[4][2], res, signe, ell]
379
}
380
381
{xpm (a, b, parms)=
382
if (parms,, parms = defaultellsyms);
383
(parms [4] * modulartofree (parms [1], parms [2], parms [3],
384
[a,b], parms [5]))[1];
385
}
386
387
388
{f(N,p)=
389
modsym(N,1); print(factor(charpoly(T(p),x)));
390
modsym(N,-1); print(factor(charpoly(T(p),x)));
391
}
392