Annotation of OpenXM_contrib/pari/src/modules/elliptic.c, Revision 1.1.1.1
1.1 maekawa 1: /********************************************************************/
2: /** **/
3: /** ELLIPTIC CURVES **/
4: /** **/
5: /********************************************************************/
6: /* $Id: elliptic.c,v 1.2 1999/09/23 17:50:57 karim Exp $ */
7: #include "pari.h"
8:
9: void
10: checkpt(GEN z)
11: {
12: if (typ(z)!=t_VEC) err(elliper1);
13: }
14:
15: long
16: checkell(GEN e)
17: {
18: long lx=lg(e);
19: if (typ(e)!=t_VEC || lx<14) err(elliper1);
20: return lx;
21: }
22:
23: void
24: checkbell(GEN e)
25: {
26: if (typ(e)!=t_VEC || lg(e)<20) err(elliper1);
27: }
28:
29: void
30: checksell(GEN e)
31: {
32: if (typ(e)!=t_VEC || lg(e)<6) err(elliper1);
33: }
34:
35: static void
36: checkch(GEN z)
37: {
38: if (typ(z)!=t_VEC || lg(z)!=5) err(elliper1);
39: }
40:
41: /* 4 X^3 + b2 X^2 + 2b4 X + b6 */
42: static GEN
43: RHSpol(GEN e)
44: {
45: GEN z = cgetg(6, t_POL); z[1] = evalsigne(1)|evallgef(6);
46: z[2] = e[8];
47: z[3] = lmul2n((GEN)e[7],1);
48: z[4] = e[6];
49: z[5] = lstoi(4); return z;
50: }
51:
52: /* x^3 + a2 x^2 + a4 x + a6 */
53: static GEN
54: ellRHS(GEN e, GEN x)
55: {
56: GEN p1;
57: p1 = gadd((GEN)e[2],x);
58: p1 = gadd((GEN)e[4], gmul(x,p1));
59: p1 = gadd((GEN)e[5], gmul(x,p1));
60: return p1;
61: }
62:
63: /* a1 x + a3 */
64: static GEN
65: ellLHS0(GEN e, GEN x)
66: {
67: return gcmp0((GEN)e[1])? (GEN)e[3]: gadd((GEN)e[3], gmul(x,(GEN)e[1]));
68: }
69:
70: static GEN
71: ellLHS0_i(GEN e, GEN x)
72: {
73: return signe(e[1])? addii((GEN)e[3], mulii(x, (GEN)e[1])): (GEN)e[3];
74: }
75:
76: /* y^2 + a1 xy + a3 y */
77: static GEN
78: ellLHS(GEN e, GEN z)
79: {
80: GEN y = (GEN)z[2];
81: return gmul(y, gadd(y, ellLHS0(e,(GEN)z[1])));
82: }
83:
84: /* 2y + a1 x + a3 */
85: static GEN
86: d_ellLHS(GEN e, GEN z)
87: {
88: return gadd(ellLHS0(e, (GEN)z[1]), gmul2n((GEN)z[2],1));
89: }
90:
91:
92: static void
93: smallinitell0(GEN x, GEN y)
94: {
95: GEN b2,b4,b6,b8,d,j,a11,a13,a33,a64,b81,b22,c4,c6;
96: long i;
97:
98: checksell(x); for (i=1; i<=5; i++) y[i]=x[i];
99:
100: b2=gadd(a11=gsqr((GEN)y[1]),gmul2n((GEN)y[2],2));
101: y[6]=(long)b2;
102:
103: b4=gadd(a13=gmul((GEN)y[1],(GEN)y[3]),gmul2n((GEN)y[4],1));
104: y[7]=(long)b4;
105:
106: b6=gadd(a33=gsqr((GEN)y[3]),a64=gmul2n((GEN)y[5],2));
107: y[8]=(long)b6;
108:
109: b81=gadd(gadd(gmul(a11,(GEN)y[5]),gmul(a64,(GEN)y[2])),gmul((GEN)y[2],a33));
110: b8=gsub(b81,gmul((GEN)y[4],gadd((GEN)y[4],a13)));
111: y[9]=(long)b8;
112:
113: c4=gadd(b22=gsqr(b2),gmulsg(-24,b4));
114: y[10]=(long)c4;
115:
116: c6=gadd(gmul(b2,gsub(gmulsg(36,b4),b22)),gmulsg(-216,b6));
117: y[11]=(long)c6;
118:
119: b81=gadd(gmul(b22,b8),gmulsg(27,gsqr(b6)));
120: d=gsub(gmul(b4,gadd(gmulsg(9,gmul(b2,b6)),gmulsg(-8,gsqr(b4)))),b81);
121: y[12]=(long)d;
122:
123: if (gcmp0(d)) err(talker,"singular curve in ellinit");
124:
125: j = gdiv(gmul(gsqr(c4),c4),d);
126: y[13]=(long)j;
127: }
128:
129: GEN
130: smallinitell(GEN x)
131: {
132: GEN y;
133: long av,tetpil;
134:
135: av=avma; y=cgetg(14,t_VEC);
136: smallinitell0(x,y); tetpil=avma;
137: return gerepile(av,tetpil,gcopy(y));
138: }
139:
140: GEN
141: ellinit0(GEN x, long flag,long prec)
142: {
143: switch(flag)
144: {
145: case 0: return initell(x,prec);
146: case 1: return smallinitell(x);
147: default: err(flagerr,"ellinit");
148: }
149: return NULL; /* not reached */
150: }
151:
152: void
153: ellprint(GEN e)
154: {
155: long av = avma;
156: long vx = fetch_var();
157: long vy = fetch_var();
158: GEN z = cgetg(3,t_VEC);
159: if (typ(e) != t_VEC || lg(e) < 6)
160: err(talker, "not an elliptic curve in ellprint");
161: z[1] = lpolx[vx]; name_var(vx, "X");
162: z[2] = lpolx[vy]; name_var(vy, "Y");
163: fprintferr("%Z = %Z\n", ellLHS(e, z), ellRHS(e, polx[vx]));
164: (void)delete_var();
165: (void)delete_var(); avma = av;
166: }
167:
168: static GEN
169: do_agm(GEN *ptx1, GEN a1, GEN b1, long prec, long sw)
170: {
171: GEN p1,r1,a,b,x,x1;
172: long G;
173:
174: x1 = gmul2n(gsub(a1,b1),-2);
175: if (gcmp0(x1))
176: err(talker,"precision too low in initell");
177: G = 6 - bit_accuracy(prec);
178: for(;;)
179: {
180: a=a1; b=b1; x=x1;
181: b1=gsqrt(gmul(a,b),prec); setsigne(b1,sw);
182: a1=gmul2n(gadd(gadd(a,b),gmul2n(b1,1)),-2);
183: r1=gsub(a1,b1);
184: p1=gsqrt(gdiv(gadd(x,r1),x),prec);
185: x1=gmul(x,gsqr(gmul2n(gaddsg(1,p1),-1)));
186: if (gexpo(r1) <= G + gexpo(b1)) break;
187: }
188: if (gprecision(x1)*2 <= (prec+2))
189: err(talker,"precision too low in initell");
190: *ptx1 = x1; return ginv(gmul2n(a1,2));
191: }
192:
193: static GEN
194: do_padic_agm(GEN *ptx1, GEN a1, GEN b1, GEN p)
195: {
196: GEN p1,r1,a,b,x,bmod1, bmod = modii((GEN)b1[4],p), x1 = *ptx1;
197:
198: if (!x1) x1 = gmul2n(gsub(a1,b1),-2);
199: for(;;)
200: {
201: a=a1; b=b1; x=x1;
202: b1=gsqrt(gmul(a,b),0); bmod1=modii((GEN)b1[4],p);
203: if (!egalii(bmod1,bmod)) b1 = gneg_i(b1);
204: a1=gmul2n(gadd(gadd(a,b),gmul2n(b1,1)),-2);
205: r1=gsub(a1,b1);
206: p1=gsqrt(gdiv(gadd(x,r1),x),0);
207: if (! gcmp1(modii((GEN)p1[4],p))) p1 = gneg_i(p1);
208: x1=gmul(x,gsqr(gmul2n(gaddsg(1,p1),-1)));
209: if (gcmp0(r1)) break;
210: }
211: *ptx1 = x1; return ginv(gmul2n(a1,2));
212: }
213:
214: static GEN
215: padic_initell(GEN y, GEN p, long prec)
216: {
217: GEN b2,b4,c4,c6,p1,p2,w,pv,a1,b1,x1,u2,q,e0,e1;
218: long i,alpha;
219:
220: if (valp(y[13]) >= 0) /* p | j */
221: err(talker,"valuation of j must be negative in p-adic ellinit");
222: if (egalii(p,gdeux))
223: err(impl,"initell for 2-adic numbers"); /* pv=stoi(4); */
224:
225: pv=p; q=ggrandocp(p,prec);
226: for (i=1; i<=5; i++) y[i]=ladd(q,(GEN)y[i]);
227: b2= (GEN)y[6];
228: b4= (GEN)y[7];
229: c4= (GEN)y[10];
230: c6= (GEN)y[11];
231: alpha=valp(c4)>>1;
232: setvalp(c4,0);
233: setvalp(c6,0); e1=gdivgs(gdiv(c6,c4),6);
234: c4=gdivgs(c4,48); c6=gdivgs(c6,864);
235: do
236: {
237: e0=e1; p2=gsqr(e0);
238: e1=gdiv(gadd(gmul2n(gmul(e0,p2),1),c6), gsub(gmulsg(3,p2),c4));
239: }
240: while (!gegal(e0,e1));
241: setvalp(e1,valp(e1)+alpha);
242:
243: e1=gsub(e1,gdivgs(b2,12));
244: w=gsqrt(gmul2n(gadd(b4,gmul(e1,gadd(b2,gmulsg(6,e1)))),1),0);
245:
246: p1=gaddgs(gdiv(gmulsg(3,e0),w),1);
247: if (valp(p1)<=0) w=gneg_i(w);
248: y[18]=(long)w;
249:
250: a1=gmul2n(gsub(w,gadd(gmulsg(3,e1),gmul2n(b2,-2))),-2);
251: b1=gmul2n(w,-1); x1=NULL;
252: u2 = do_padic_agm(&x1,a1,b1,pv);
253:
254: w=gaddsg(1,ginv(gmul2n(gmul(u2,x1),1)));
255: q=ginv(gadd(w,gsqrt(gaddgs(gsqr(w),-1),0)));
256: if (valp(q)<0) q=ginv(q);
257:
258: p1=cgetg(2,t_VEC); p1[1]=(long)e1;
259: y[14]=(long)p1;
260: y[15]=(long)u2;
261: y[16] = (kronecker((GEN)u2[4],p) <= 0 || (valp(u2)&1))? zero: lsqrt(u2,0);
262: y[17]=(long)q;
263: y[19]=zero; return y;
264: }
265:
266: static int
267: invcmp(GEN x, GEN y) { return -gcmp(x,y); }
268:
269: static GEN
270: initell0(GEN x, long prec)
271: {
272: GEN y,b2,b4,D,p1,p2,p,u,w,a1,b1,x1,u2,q,e1,pi,pi2,tau,w2;
273: long ty,i,e,sw;
274:
275: y=cgetg(20,t_VEC);
276: smallinitell0(x,y);
277:
278: e = BIGINT; p = NULL;
279: for (i=1; i<=5; i++)
280: {
281: q = (GEN)y[i];
282: if (typ(q)==t_PADIC)
283: {
284: long e2 = signe(q[4])? precp(q)+valp(q): valp(q);
285: if (e2 < e) e = e2;
286: if (!p) p = (GEN)q[2];
287: else if (!egalii(p,(GEN)q[2]))
288: err(talker,"incompatible p-adic numbers in initell");
289: }
290: }
291: if (e<BIGINT) return padic_initell(y,p,e);
292:
293: b2= (GEN)y[6];
294: b4= (GEN)y[7];
295: D = (GEN)y[12]; ty = typ(D);
296: if (!prec || !is_const_t(ty) || ty==t_INTMOD)
297: { y[14]=y[15]=y[16]=y[17]=y[18]=y[19]=zero; return y; }
298:
299: pi=mppi(prec); pi2=gmul2n(pi,1);
300: p1 = roots(RHSpol(y),prec);
301: if (gsigne(D) < 0) p1[1] = lreal((GEN)p1[1]);
302: else /* sort roots in decreasing order */
303: p1 = gen_sort(greal(p1), 0, invcmp);
304: y[14]=(long)p1;
305:
306: e1=(GEN)p1[1];
307: w = gsqrt(gmul2n(gadd(b4,gmul(e1,gadd(b2,gmulsg(6,e1)))),1),prec);
308: p2=gadd(gmulsg(3,e1), gmul2n(b2,-2));
309: if (gsigne(p2)>0) w = gneg_i(w);
310: a1=gmul2n(gsub(w,p2),-2);
311: b1=gmul2n(w,-1); sw=signe(w);
312: u2 = do_agm(&x1,a1,b1,prec,sw);
313:
314: w = gaddsg(1,ginv(gmul2n(gmul(u2,x1),1)));
315: q = gsqrt(gaddgs(gsqr(w),-1),prec);
316: if (gsigne(greal(w))>0)
317: q = ginv(gadd(w,q));
318: else
319: q = gsub(w,q);
320: if (gexpo(q)>=0) q=ginv(q);
321: tau=gmul(gdiv(glog(q,prec),pi2), gneg_i(gi));
322:
323: y[19]=lmul(gmul(gsqr(pi2),gabs(u2,prec)),gimag(tau));
324: u=gmul(pi2,gsqrt(gneg_i(u2),prec)); w2=gmul(tau,u);
325: if (sw<0)
326: {
327: y[15]=(long)u;
328: q=gsqrt(q,prec);
329: }
330: else
331: {
332: y[15]=lmul2n(gabs((GEN)w2[1],prec),1);
333: q=gexp(gmul2n(gmul(gmul(pi2,gi),gdiv(w2,(GEN)y[15])), -1), prec);
334: }
335: y[16]=(long)w2;
336: p1 = gdiv(gsqr(pi),gmulsg(6,(GEN)y[15]));
337: p2 = thetanullk(q,1,prec);
338: if (gcmp0(p2)) err(talker,"precision too low in initell");
339: y[17]=lmul(p1,gdiv(thetanullk(q,3,prec),p2));
340: y[18]=ldiv(gsub(gmul((GEN)y[17],(GEN)y[16]),gmul(gi,pi)),(GEN)y[15]);
341: return y;
342: }
343:
344: GEN
345: initell(GEN x, long prec)
346: {
347: long av=avma;
348: return gerepileupto(av, gcopy(initell0(x,prec)));
349: }
350:
351: GEN
352: coordch(GEN e, GEN ch)
353: {
354: GEN y,p1,p2,v,v2,v3,v4,v6,r,s,t,u;
355: long av,tetpil,i,lx = checkell(e);
356:
357: checkch(ch);
358: u=(GEN)ch[1]; r=(GEN)ch[2]; s=(GEN)ch[3]; t=(GEN)ch[4];
359: av=avma; y=cgetg(lx,t_VEC);
360: v=ginv(u); v2=gsqr(v); v3=gmul(v,v2);v4=gsqr(v2); v6=gsqr(v3);
361: y[1] = lmul(v,gadd((GEN)e[1],gmul2n(s,1)));
362: y[2] = lmul(v2,gsub(gadd((GEN)e[2],gmulsg(3,r)),gmul(s,gadd((GEN)e[1],s))));
363: p2 = ellLHS0(e,r);
364: p1 = gadd(gmul2n(t,1), p2);
365: y[3] = lmul(v3,p1);
366: p1 = gsub((GEN)e[4],gadd(gmul(t,(GEN)e[1]),gmul(s,p1)));
367: y[4] = lmul(v4,gadd(p1,gmul(r,gadd(gmul2n((GEN)e[2],1),gmulsg(3,r)))));
368: p2 = gmul(t,gadd(t, p2));
369: y[5] = lmul(v6,gsub(ellRHS(e,r),p2));
370: y[6] = lmul(v2,gadd((GEN)e[6],gmulsg(12,r)));
371: y[7] = lmul(v4,gadd((GEN)e[7],gmul(r,gadd((GEN)e[6],gmulsg(6,r)))));
372: y[8] = lmul(v6,gadd((GEN)e[8],gmul(r,gadd(gmul2n((GEN)e[7],1),gmul(r,gadd((GEN)e[6],gmul2n(r,2)))))));
373: p1 = gadd(gmulsg(3,(GEN)e[7]),gmul(r,gadd((GEN)e[6],gmulsg(3,r))));
374: y[9] = lmul(gsqr(v4),gadd((GEN)e[9],gmul(r,gadd(gmulsg(3,(GEN)e[8]),gmul(r,p1)))));
375: y[10] = lmul(v4,(GEN)e[10]);
376: y[11] = lmul(v6,(GEN)e[11]);
377: y[12] = lmul(gsqr(v6),(GEN)e[12]);
378: y[13] = e[13];
379: if (lx>14)
380: {
381: p1=(GEN)e[14];
382: if (gcmp0(p1))
383: {
384: y[14] = y[15] = y[16] = y[17] = y[18] = y[19] = zero;
385: }
386: else
387: {
388: if (typ(e[1])==t_PADIC)
389: {
390: p2=cgetg(2,t_VEC); p2[1]=lmul(v2,gsub((GEN)p1[1],r));
391: y[14]=(long)p2;
392: y[15]=lmul(gsqr(u),(GEN)e[15]);
393: y[16]=lmul(u,(GEN)e[16]);
394: /* FIXME: how do q and w change ??? */
395: y[17]=e[17];
396: y[18]=e[18];
397: y[19]=zero;
398: }
399: else
400: {
401: p2=cgetg(4,t_COL);
402: for (i=1; i<=3; i++) p2[i]=lmul(v2,gsub((GEN)p1[i],r));
403: y[14]=(long)p2;
404: y[15]=lmul(u,(GEN)e[15]);
405: y[16]=lmul(u,(GEN)e[16]);
406: y[17]=ldiv((GEN)e[17],u);
407: y[18]=ldiv((GEN)e[18],u);
408: y[19]=lmul(gsqr(u),(GEN)e[19]);
409: }
410: }
411: }
412: tetpil=avma; return gerepile(av,tetpil,gcopy(y));
413: }
414:
415: static GEN
416: pointch0(GEN x, GEN v2, GEN v3, GEN mor, GEN s, GEN t)
417: {
418: GEN p1,z;
419:
420: if (lg(x) < 3) return x;
421:
422: z = cgetg(3,t_VEC); p1=gadd((GEN)x[1],mor);
423: z[1] = lmul(v2,p1);
424: z[2] = lmul(v3,gsub((GEN)x[2],gadd(gmul(s,p1),t)));
425: return z;
426: }
427:
428: GEN
429: pointch(GEN x, GEN ch)
430: {
431: GEN y,v,v2,v3,mor,r,s,t,u;
432: long av,tetpil,tx,lx=lg(x),i;
433:
434: checkpt(x); checkch(ch);
435: if (lx < 2) return gcopy(x);
436: av=avma; u=(GEN)ch[1]; r=(GEN)ch[2]; s=(GEN)ch[3]; t=(GEN)ch[4];
437: tx=typ(x[1]); v=ginv(u); v2=gsqr(v); v3=gmul(v,v2); mor=gneg_i(r);
438: if (is_matvec_t(tx))
439: {
440: y=cgetg(lx,tx);
441: for (i=1; i<lx; i++)
442: y[i]=(long) pointch0((GEN)x[i],v2,v3,mor,s,t);
443: }
444: else
445: y = pointch0(x,v2,v3,mor,s,t);
446: tetpil=avma; return gerepile(av,tetpil,gcopy(y));
447: }
448:
449: static long
450: ellexpo(GEN e)
451: {
452: long i,k2, k = gexpo((GEN)e[1]);
453: for (i=2; i<6; i++)
454: {
455: k2 = gexpo((GEN)e[i]);
456: if (k<k2) k = k2;
457: }
458: return k;
459: }
460:
461: /* Exactness of lhs and rhs in the following depends in non-obvious ways
462: on the coeffs of the curve as well as on the components of the point z.
463: Thus if e is exact, with a1==0, and z has exact y coordinate only, the
464: lhs will be exact but the rhs won't. */
465: int
466: oncurve(GEN e, GEN z)
467: {
468: GEN p1,p2,x;
469: long av=avma,p,q;
470:
471: checksell(e); checkpt(z); if (lg(z)<3) return 1; /* oo */
472: p1 = ellLHS(e,z);
473: p2 = ellRHS(e,(GEN)z[1]); x = gsub(p1,p2);
474: if (gcmp0(x)) { avma=av; return 1; }
475: p = precision(p1);
476: q = precision(p2);
477: if (!p && !q) { avma=av; return 0; } /* both of p1, p2 are exact */
478: if (!q || (p && p < q)) q = p; /* min among nonzero elts of {p,q} */
479: /* the constant 0.93 is arbitrary */
480: q = (gexpo(x)-ellexpo(e) < - bit_accuracy(q) * 0.93);
481: avma = av; return q;
482: }
483:
484: GEN
485: addell(GEN e, GEN z1, GEN z2)
486: {
487: GEN p1,p2,x,y,x1,x2,y1,y2;
488: long av=avma,tetpil;
489:
490: checksell(e); checkpt(z1); checkpt(z2);
491: if (lg(z1)<3) return gcopy(z2);
492: if (lg(z2)<3) return gcopy(z1);
493:
494: x1=(GEN)z1[1]; y1=(GEN)z1[2];
495: x2=(GEN)z2[1]; y2=(GEN)z2[2];
496: if (x1 == x2 || gegal(x1,x2))
497: { /* y1 = y2 or -LHS0-y2 */
498: if (y1 != y2)
499: {
500: int eq;
501: if (precision(y1) || precision(y2))
502: eq = (gexpo(gadd(ellLHS0(e,x1),gadd(y1,y2))) >= gexpo(y1));
503: else
504: eq = gegal(y1,y2);
505: if (!eq) { avma=av; y=cgetg(2,t_VEC); y[1]=zero; return y; }
506: }
507: p2 = d_ellLHS(e,z1);
508: if (gcmp0(p2)) { avma=av; y=cgetg(2,t_VEC); y[1]=zero; return y; }
509: p1 = gadd(gsub((GEN)e[4],gmul((GEN)e[1],y1)),
510: gmul(x1,gadd(gmul2n((GEN)e[2],1),gmulsg(3,x1))));
511: }
512: else { p1=gsub(y2,y1); p2=gsub(x2,x1); }
513: p1 = gdiv(p1,p2);
514: x = gsub(gmul(p1,gadd(p1,(GEN)e[1])), gadd(gadd(x1,x2),(GEN)e[2]));
515: y = gadd(gadd(y1, ellLHS0(e,x)), gmul(p1,gsub(x,x1)));
516: tetpil=avma; p1=cgetg(3,t_VEC); p1[1]=lcopy(x); p1[2]=lneg(y);
517: return gerepile(av,tetpil,p1);
518: }
519:
520: static GEN
521: invell(GEN e, GEN z)
522: {
523: GEN p1;
524:
525: if (lg(z)<3) return z;
526: p1=cgetg(3,t_VEC); p1[1]=z[1];
527: p1[2]=(long)gneg_i(gadd((GEN)z[2], ellLHS0(e,(GEN)z[1])));
528: return p1;
529: }
530:
531: GEN
532: subell(GEN e, GEN z1, GEN z2)
533: {
534: long av=avma,tetpil;
535:
536: checksell(e); checkpt(z2);
537: z2=invell(e,z2); tetpil=avma;
538: return gerepile(av,tetpil,addell(e,z1,z2));
539: }
540:
541: GEN
542: ordell(GEN e, GEN x, long prec)
543: {
544: long av=avma,td,i,lx,tx=typ(x);
545: GEN D,a,b,d,pd,p1,y;
546:
547: checksell(e);
548: if (is_matvec_t(tx))
549: {
550: lx=lg(x); y=cgetg(lx,tx);
551: for (i=1; i<lx; i++) y[i]=(long)ordell(e,(GEN)x[i],prec);
552: return y;
553: }
554:
555: a=ellRHS(e,x);
556: b=ellLHS0(e,x); /* y*(y+b) = a */
557: D=gadd(gsqr(b),gmul2n(a,2)); td=typ(D);
558: if (gcmp0(D))
559: {
560: b = gneg_i(b);
561: y = cgetg(2,t_VEC);
562: if (td == t_INTMOD && egalii((GEN)D[1], gdeux))
563: y[1] = (long)gmodulss(gcmp0(a)?0:1, 2);
564: else
565: y[1] = lmul2n(b,-1);
566: return gerepileupto(av,y);
567: }
568:
569: if (td==t_INT || is_frac_t(td))
570: {
571: pd = (td==t_INT)? NULL: (GEN)D[2];
572: if (pd) D = mulii((GEN)D[1],pd);
573: if (!carrecomplet(D,&d)) { avma=av; return cgetg(1,t_VEC); }
574: if (pd) d = gdiv(d,pd);
575: }
576: else
577: {
578: if (td==t_INTMOD)
579: {
580: if (egalii((GEN)D[1],gdeux))
581: {
582: avma=av;
583: if (!gcmp0(a)) return cgetg(1,t_VEC);
584: y = cgetg(3,t_VEC);
585: y[1] = (long)gmodulss(0,2);
586: y[2] = (long)gmodulss(1,2); return y;
587: }
588: if (kronecker((GEN)D[2],(GEN)D[1]) == -1)
589: { avma=av; return cgetg(1,t_VEC); }
590: }
591: d = gsqrt(D,prec);
592: }
593: p1=gsub(d,b); y = cgetg(3,t_VEC);
594: y[1] = lmul2n(p1,-1);
595: y[2] = lsub((GEN)y[1],d);
596: return gerepileupto(av,y);
597: }
598:
599: static GEN
600: CM_powell(GEN e, GEN z, GEN n)
601: {
602: GEN x,y,p0,p1,q0,q1,p2,q2,z1,z2,pol,grdx;
603: long av=avma,tetpil,ln,ep,vn;
604:
605: if (lg(z)<3) return gcopy(z);
606: pol=(GEN)n[1];
607: if (signe(discsr(pol))>=0)
608: err(talker,"not a negative quadratic discriminant in CM");
609: if (!gcmp1(denom((GEN)n[2])) || !gcmp1(denom((GEN)n[3])))
610: err(impl,"powell for nonintegral CM exponent");
611:
612: p1=gaddgs(gmul2n(gnorm(n),2),4);
613: if (gcmpgs(p1,(((ulong)MAXULONG)>>1)) > 0)
614: err(talker,"norm too large in CM");
615: ln=itos(p1); vn=(ln-4)>>2;
616: z1 = weipell(e,ln);
617: z2 = gsubst(z1,0,gmul(n,polx[0]));
618: grdx=gadd((GEN)z[1],gdivgs((GEN)e[6],12));
619: p0=gzero; p1=gun;
620: q0=gun; q1=gzero;
621: do
622: {
623: GEN ss=gzero;
624: do
625: {
626: ep=(-valp(z2))>>1; ss=gadd(ss,gmul((GEN)z2[2],gpuigs(polx[0],ep)));
627: z2=gsub(z2,gmul((GEN)z2[2],gpuigs(z1,ep)));
628: }
629: while (valp(z2)<=0);
630: p2=gadd(p0,gmul(ss,p1)); p0=p1; p1=p2;
631: q2=gadd(q0,gmul(ss,q1)); q0=q1; q1=q2;
632: if (!signe(z2)) break;
633: z2=ginv(z2);
634: }
635: while (lgef(p1)-3 < vn);
636: if (lgef(p1)-3 > vn || signe(z2))
637: err(talker,"not a complex multiplication in powell");
638: x=gdiv(p1,q1); y=gdiv(deriv(x,0),n);
639: x=gsub(poleval(x,grdx), gdivgs((GEN)e[6],12));
640: y=gsub(gmul(d_ellLHS(e,z),poleval(y,grdx)), ellLHS0(e,x));
641: tetpil=avma; z=cgetg(3,t_VEC); z[1]=lcopy(x); z[2]=lmul2n(y,-1);
642: return gerepile(av,tetpil,z);
643: }
644:
645: GEN
646: powell(GEN e, GEN z, GEN n)
647: {
648: GEN y;
649: long av=avma,i,j,tetpil,s;
650: ulong m;
651:
652: checksell(e); checkpt(z);
653: if (typ(n)==t_QUAD) return CM_powell(e,z,n);
654: if (typ(n)!=t_INT)
655: err(impl,"powell for nonintegral or non CM exponents");
656: if (lg(z)<3) return gcopy(z);
657: s=signe(n);
658: if (!s) { y=cgetg(2,t_VEC); y[1]=zero; return y; }
659: if (s < 0) { n=negi(n); z = invell(e,z); }
660: if (is_pm1(n)) { tetpil=avma; return gerepile(av,tetpil,gcopy(z)); }
661:
662: y=cgetg(2,t_VEC); y[1]=zero;
663: for (i=lgefint(n)-1; i>2; i--)
664: for (m=n[i],j=0; j<BITS_IN_LONG; j++,m>>=1)
665: {
666: if (m&1) y = addell(e,y,z);
667: z = addell(e,z,z);
668: }
669: for (m=n[2]; m>1; m>>=1)
670: {
671: if (m&1) y = addell(e,y,z);
672: z = addell(e,z,z);
673: }
674: tetpil=avma; return gerepile(av,tetpil,addell(e,y,z));
675: }
676:
677: GEN
678: mathell(GEN e, GEN x, long prec)
679: {
680: GEN y,p1,p2, *pdiag;
681: long av=avma,tetpil,lx=lg(x),i,j,tx=typ(x);
682:
683: if (!is_vec_t(tx)) err(elliper1);
684: lx=lg(x); y=cgetg(lx,t_MAT); pdiag=(GEN*) new_chunk(lx);
685: for (i=1; i<lx; i++)
686: {
687: pdiag[i]=ghell(e,(GEN)x[i],prec);
688: y[i]=lgetg(lx,t_COL);
689: }
690: for (i=1; i<lx; i++)
691: {
692: p1=(GEN)y[i]; p1[i]=lmul2n(pdiag[i],1);
693: for (j=i+1; j<lx; j++)
694: {
695: p2=ghell(e,addell(e,(GEN)x[i],(GEN)x[j]),prec);
696: p2=gsub(p2, gadd(pdiag[i],pdiag[j]));
697: p1[j]=(long)p2; coeff(y,i,j)=(long)p2;
698: }
699: }
700: tetpil=avma; return gerepile(av,tetpil,gcopy(y));
701: }
702:
703: static GEN
704: bilhells(GEN e, GEN z1, GEN z2, GEN h2, long prec)
705: {
706: long lz1=lg(z1),tx,av=avma,tetpil,i;
707: GEN y,p1,p2;
708:
709: if (lz1==1) return cgetg(1,typ(z1));
710:
711: tx=typ(z1[1]);
712: if (!is_matvec_t(tx))
713: {
714: p1 = ghell(e,addell(e,z1,z2),prec);
715: p2 = gadd(ghell(e,z1,prec),h2);
716: tetpil=avma; return gerepile(av,tetpil,gsub(p1,p2));
717: }
718: y=cgetg(lz1,typ(z1));
719: for (i=1; i<lz1; i++)
720: y[i]=(long)bilhells(e,(GEN)z1[i],z2,h2,prec);
721: return y;
722: }
723:
724: GEN
725: bilhell(GEN e, GEN z1, GEN z2, long prec)
726: {
727: GEN p1,h2;
728: long av=avma,tetpil,tz1=typ(z1),tz2=typ(z2);
729:
730: if (!is_matvec_t(tz1) || !is_matvec_t(tz2)) err(elliper1);
731: if (lg(z1)==1) return cgetg(1,tz1);
732: if (lg(z2)==1) return cgetg(1,tz2);
733:
734: tz1=typ(z1[1]); tz2=typ(z2[1]);
735: if (is_matvec_t(tz2))
736: {
737: if (is_matvec_t(tz1))
738: err(talker,"two vector/matrix types in bilhell");
739: p1=z1; z1=z2; z2=p1;
740: }
741: h2=ghell(e,z2,prec); tetpil=avma;
742: return gerepile(av,tetpil,bilhells(e,z1,z2,h2,prec));
743: }
744:
745: static GEN
746: new_coords(GEN e, GEN x, GEN *pta, GEN *ptb, long prec)
747: {
748: GEN a,b,r0,r1,p1,p2,w, e1 = gmael(e,14,1), b2 = (GEN)e[6];
749: long ty = typ(e[12]);
750:
751: r0 = gmul2n(b2,-2);
752: p2 = gadd(gmulsg(3,e1),r0);
753: if (ty == t_PADIC)
754: w = (GEN)e[18];
755: else
756: {
757: GEN b4 = (GEN)e[7];
758: if (!is_const_t(ty)) err(typeer,"zell");
759:
760: w = gsqrt(gmul2n(gadd(b4, gmul(e1,gadd(b2,gmulsg(6,e1)))),1),prec);
761: if (gsigne(greal(p2)) > 0) w = gneg_i(w);
762: }
763: a = gmul2n(gsub(w,p2),-2);
764: b = gmul2n(w,-1);
765: r1 = gsub(a,b);
766: p1 = gadd(x, gmul2n(gadd(e1,r0),-1));
767: if (gcmp0(p1)) p1=gsqrt(gneg_i(gmul(a,r1)),prec);
768: else
769: {
770: GEN delta=gdiv(gmul(a,r1),gsqr(p1));
771: p1=gmul2n(gmul(p1,gaddsg(1,gsqrt(gsubsg(1,gmul2n(delta,2)),prec))),-1);
772: }
773: *pta = a; *ptb = b;
774: return gmul(p1,gsqr(gmul2n(gaddsg(1,gsqrt(gdiv(gadd(p1,r1),p1),prec)),-1)));
775: }
776:
777: GEN
778: zell(GEN e, GEN z, long prec)
779: {
780: long av=avma,ty,sw,fl;
781: GEN t,u,p1,p2,r1,a,b,x1,u2,D = (GEN)e[12];
782:
783: checkbell(e);
784: if (!oncurve(e,z)) err(heller1);
785: ty=typ(D);
786: if (ty==t_INTMOD) err(typeer,"zell");
787: if (lg(z)<3) return (ty==t_PADIC)? gun: gzero;
788:
789: x1 = new_coords(e,(GEN)z[1],&a,&b,prec);
790: if (ty==t_PADIC)
791: {
792: u2 = do_padic_agm(&x1,a,b,(GEN)D[2]);
793: if (!gcmp0((GEN)e[16]))
794: {
795: t=gsqrt(gaddsg(1,gdiv(x1,a)),prec);
796: t=gdiv(gaddsg(-1,t),gaddsg(1,t));
797: }
798: else t=gaddsg(2,ginv(gmul(u2,x1)));
799: return gerepileupto(av,t);
800: }
801:
802: sw=gsigne(greal(b)); fl=0;
803: for(;;) /* agm */
804: {
805: GEN a0=a, b0=b, x0=x1;
806:
807: b=gsqrt(gmul(a0,b0),prec);
808: if (gsigne(greal(b)) != sw) b = gneg_i(b);
809: a=gmul2n(gadd(gadd(a0,b0),gmul2n(b,1)),-2);
810: r1=gsub(a,b);
811: p1=gsqrt(gdiv(gadd(x0,r1),x0),prec);
812: x1=gmul(x0,gsqr(gmul2n(gaddsg(1,p1),-1)));
813: if (gexpo(gsub(x1,x0)) < gexpo(x1) - bit_accuracy(prec) + 5)
814: {
815: if (fl) break;
816: fl = 1;
817: }
818: else fl = 0;
819: }
820: u=gdiv(x1,a); t=gaddsg(1,u);
821: if (gexpo(t) < 5 - bit_accuracy(prec))
822: t = negi(gun);
823: else
824: t = gdiv(u,gsqr(gaddsg(1,gsqrt(t,prec))));
825: u=gsqrt(ginv(gmul2n(a,2)),prec);
826: t=glog(t,prec); t=gmul(t,u);
827:
828: /* which square root? test the reciprocal function (pointell) */
829: if (!gcmp0(t))
830: {
831: GEN x1;
832: long bad;
833:
834: u = pointell(e,t,3); /* we don't need much precision */
835: /* Either z = u (ok: keep t), or z = invell(e,u) (bad: t <-- -t) */
836: x1 = gsub(z,u); bad = (gexpo((GEN)x1[1]) >= gexpo((GEN)u[1])
837: || gexpo((GEN)x1[2]) >= gexpo((GEN)u[2]));
838: if (bad) t = gneg(t);
839: if (DEBUGLEVEL)
840: {
841: if (DEBUGLEVEL>4)
842: {
843: fprintferr(" z = %Z\n",z);
844: fprintferr(" u = %Z\n",u);
845: fprintferr(" x1 = %Z\n",x1);
846: }
847: fprintferr("ellpointtoz: %s square root\n", bad? "bad": "good");
848: flusherr();
849: }
850: }
851: /* send t to the fundamental domain if necessary */
852: p2 = gdiv(gimag(t),gmael(e,16,2));
853: p1 = gsub(p2, gmul2n(gun,-2));
854: if (gcmp(gabs(p1,prec),ghalf) >= 0)
855: t = gsub(t, gmul((GEN)e[16],gfloor(gadd(p2,dbltor(0.1)))));
856: if (gsigne(greal(t)) < 0) t = gadd(t,(GEN)e[15]);
857: return gerepileupto(av,t);
858: }
859:
860: /* compute gamma in SL_2(Z) and t'=gamma(t) so that t' is in the usual
861: fundamental domain. Internal function no check, no garbage. */
862: static GEN
863: getgamma(GEN t)
864: {
865: GEN a,b,c,d,n,m,y,p1,unapprox;
866:
867: unapprox=gsub(gun,gpuigs(stoi(10),-8));
868: a=d=gun;b=c=gzero;
869: for(;;)
870: {
871: n=ground(greal(t));
872: if (signe(n))
873: {
874: n=negi(n); t=gadd(t,n);
875: a=addii(a,mulii(n,c));
876: b=addii(b,mulii(n,d));
877: }
878: m=gnorm(t); if (gcmp(m,unapprox)>=0) break;
879: t=gneg_i(gdiv(gconj(t),m));
880: p1=negi(c); c=a; a=p1;
881: p1=negi(d); d=b; b=p1;
882: }
883: y=cgetg(3,t_VEC);
884: m=cgetg(3,t_MAT); y[1]=(long)m;
885: p1=cgetg(3,t_COL); m[1]=(long)p1;
886: p1[1]=(long)a; p1[2]=(long)c;
887: p1=cgetg(3,t_COL); m[2]=(long)p1;
888: p1[1]=(long)b; p1[2]=(long)d;
889: y[2]=(long)t;
890: return y;
891: }
892:
893: static int
894: get_periods(GEN e, GEN *om1, GEN *om2)
895: {
896: long tx = typ(e);
897: if (is_vec_t(tx))
898: switch(lg(e))
899: {
900: case 3: *om1=(GEN)e[1]; *om2=(GEN)e[2]; return 1;
901: case 20: *om1=(GEN)e[16]; *om2=(GEN)e[15]; return 1;
902: }
903: return 0;
904: }
905:
906: /* computes the numerical values of eisenstein series. k is equal to a positive
907: even integer. If k=4 or 6, computes g2 or g3. If k=2, or k>6 even,
908: compute (2iPi/om2)^k*(1+2/zeta(1-k)*sum(n>=1,n^(k-1)q^n/(1-q^n)) with no
909: constant factor. */
910: GEN
911: elleisnum(GEN om, long k, long flag, long prec)
912: {
913: long av=avma,lim,av1,fl,si;
914: GEN om1,om2,p1,pii2,tau,q,y,qn,v,ga,court,asub;
915:
916: if (k%2 || k<=0) err(talker,"k not a positive even integer in elleisnum");
917: if (!get_periods(om, &om1, &om2)) err(typeer,"elleisnum");
918: p1=mppi(prec); setexpo(p1,2);
919: pii2=cgetg(3,t_COMPLEX); pii2[1]=zero; pii2[2]=(long)p1;
920: tau=gdiv(om1,om2); si=gsigne(gimag(tau));
921: if (si==0)
922: err(talker,"omega1 and omega2 are R-linearly dependent in elleisnum");
923: if (si<0) { p1=om1; om1=om2; om2=p1; tau=ginv(tau); }
924: v=getgamma(tau); tau=(GEN)v[2]; ga=(GEN)v[1];
925: if (k==2) asub=gdiv(gmul(pii2,gmulsg(12,gcoeff(ga,2,1))),om2);
926: om2=gadd(gmul(gcoeff(ga,2,1),om1),gmul(gcoeff(ga,2,2),om2));
927: if (k==2) asub=gdiv(asub,om2);
928: q=gexp(gmul(pii2,tau),prec);
929: y=gzero; court=stoi(3);
930: av1=avma; lim=stack_lim(av1,1); qn=gun; court[2]=0;
931: do
932: {
933: court[2]++; qn=gmul(q,qn);
934: p1=gdiv(gmul(gpuigs(court,k-1),qn),gsub(gun,qn));
935: y=gadd(y,p1);
936: fl=(gexpo(p1) > - bit_accuracy(prec) - 5);
937: if (low_stack(lim, stack_lim(av1,1)))
938: {
939: GEN *gptr[2]; gptr[0]=&y; gptr[1]=&qn;
940: if(DEBUGMEM>1) err(warnmem,"elleisnum");
941: gerepilemany(av1,gptr,2);
942: }
943: }
944: while (fl);
945:
946: y=gadd(gun,gmul(gdiv(gdeux,gzeta(stoi(1-k),prec)),y));
947: p1=gpuigs(gdiv(pii2,om2),k);
948: y = gmul(p1,y);
949: if (k==2) y=gsub(y,asub);
950: else if (k==4 && flag) y=gdivgs(y,12);
951: else if (k==6 && flag) y=gdivgs(y,216);
952: return gerepileupto(av,y);
953: }
954:
955: /* compute eta1, eta2 */
956:
957: GEN
958: elleta(GEN om, long prec)
959: {
960: long av=avma,tetpil;
961: GEN e2,p1,pii2,y1,y2,y;
962:
963: e2=gdivgs(elleisnum(om,2,0,prec),12);
964: p1=mppi(prec); setexpo(p1,2);
965: pii2=cgetg(3,t_COMPLEX); pii2[1]=zero; pii2[2]=(long)p1;
966: y2=gmul((GEN)om[2],e2);
967: y1=gadd(gdiv(pii2,(GEN)om[2]),gmul((GEN)om[1],e2));
968: tetpil=avma;
969: y=cgetg(3,t_VEC); y[1]=lneg(y1); y[2]=lneg(y2);
970: return gerepile(av,tetpil,y);
971: }
972:
973: /* computes the numerical value of wp(z | om1 Z + om2 Z),
974: If flall=1, compute also wp'. Reduce to the fundamental domain first. */
975: static GEN
976: weipellnumall(GEN om1, GEN om2, GEN z, long flall, long prec)
977: {
978: long av=avma,tetpil,lim,av1,si,toadd;
979: GEN p1,pii2,pii4,a,tau,q,u,y,yp,u1,u2,qn,v,ga;
980:
981: p1=mppi(prec); setexpo(p1,2);
982: pii2=cgetg(3,t_COMPLEX); pii2[1]=zero; pii2[2]=(long)p1;
983: tau=gdiv(om1,om2); si=gsigne(gimag(tau));
984: if (si==0)
985: err(talker,"omega1 and omega2 are R-linearly dependent in ellwpnum");
986: if (si<0) { p1=om1; om1=om2; om2=p1; tau=ginv(tau); }
987: v=getgamma(tau); tau=(GEN)v[2]; ga=(GEN)v[1];
988: om2=gadd(gmul(gcoeff(ga,2,1),om1),gmul(gcoeff(ga,2,2),om2));
989: z=gdiv(z,om2);
990: a=ground(gdiv(gimag(z),gimag(tau))); z=gsub(z,gmul(a,tau));
991: a=ground(greal(z)); z=gsub(z,a);
992: if (gexpo(z) < 5 - bit_accuracy(prec))
993: {
994: avma=av; v=cgetg(2,t_VEC); v[1]=zero; return v;
995: }
996:
997: q=gexp(gmul(pii2,tau),prec);
998: u=gexp(gmul(pii2,z),prec);
999: u1=gsub(gun,u); u2=gsqr(u1);
1000: y=gadd(gdivgs(gun,12),gdiv(u,u2));
1001: if (flall) yp=gdiv(gadd(gun,u),gmul(u1,u2));
1002: toadd=(long)ceil(9.065*gtodouble(gimag(z)));
1003:
1004: av1=avma; lim=stack_lim(av1,1); qn=q;
1005: do
1006: {
1007: GEN p2,qnu,qnu1,qnu2,qnu3,qnu4;
1008:
1009: qnu=gmul(qn,u); qnu1=gsub(gun,qnu); qnu2=gsqr(qnu1);
1010: qnu3=gsub(qn,u); qnu4=gsqr(qnu3);
1011: p1=gsub(gmul(u,gadd(ginv(qnu2),ginv(qnu4))),
1012: gmul2n(ginv(gsqr(gsub(gun,qn))),1));
1013: p1=gmul(qn,p1);
1014: y=gadd(y,p1);
1015: if (flall)
1016: {
1017: p2=gadd(gdiv(gadd(gun,qnu),gmul(qnu1,qnu2)),
1018: gdiv(gadd(qn,u),gmul(qnu3,qnu4)));
1019: p2=gmul(qn,p2);
1020: yp=gadd(yp,p2);
1021: }
1022: qn=gmul(q,qn);
1023: if (low_stack(lim, stack_lim(av1,1)))
1024: {
1025: GEN *gptr[3]; gptr[0]=&y; gptr[1]=&qn; gptr[2]=&yp;
1026: if(DEBUGMEM>1) err(warnmem,"weipellnum");
1027: gerepilemany(av1,gptr,flall?3:2);
1028: }
1029: }
1030: while (gexpo(qn) > - bit_accuracy(prec) - 5 - toadd);
1031:
1032: pii2=gdiv(pii2,om2);
1033: pii4=gsqr(pii2);
1034: y = gmul(pii4,y);
1035: if (flall) yp=gmul(u,gmul(gmul(pii4,pii2),yp));
1036: tetpil=avma;
1037: if (flall) { v=cgetg(3,t_VEC); v[1]=lcopy(y); v[2]=lmul2n(yp,-1); }
1038: else v=gcopy(y);
1039: return gerepile(av,tetpil,v);
1040: }
1041:
1042: GEN
1043: ellzeta(GEN om, GEN z, long prec)
1044: {
1045: long av=avma,tetpil,lim,av1,si,toadd;
1046: GEN zinit,om1,om2,p1,pii2,tau,q,u,y,u1,qn,v,ga,x1,x2,et;
1047:
1048: if (!get_periods(om, &om1, &om2)) err(typeer,"ellzeta");
1049: p1=mppi(prec); setexpo(p1,2);
1050: pii2=cgetg(3,t_COMPLEX); pii2[1]=zero; pii2[2]=(long)p1;
1051: tau=gdiv(om1,om2); si=gsigne(gimag(tau));
1052: if (si==0)
1053: err(talker,"omega1 and omega2 are R-linearly dependent in ellzeta");
1054: if (si<0) { p1=om1; om1=om2; om2=p1; tau=ginv(tau); }
1055: v=getgamma(tau); tau=(GEN)v[2]; ga=(GEN)v[1];
1056: om2=gadd(gmul(gcoeff(ga,2,1),om1),gmul(gcoeff(ga,2,2),om2));
1057: om1=gmul(tau,om2); om=cgetg(3,t_VEC); om[1]=(long)om1; om[2]=(long)om2;
1058: z=gdiv(z,om2);
1059:
1060: x1=ground(gdiv(gimag(z),gimag(tau))); z=gsub(z,gmul(x1,tau));
1061: x2=ground(greal(z)); z=gsub(z,x2); zinit=gmul(z,om2);
1062: et=elleta(om,prec);
1063: et=gadd(gmul(x1,(GEN)et[1]),gmul(x2,(GEN)et[2]));
1064: if (gexpo(z) < 5 - bit_accuracy(prec))
1065: {
1066: p1=ginv(zinit); tetpil=avma; return gerepile(av,tetpil,gadd(p1,et));
1067: }
1068: q=gexp(gmul(pii2,tau),prec);
1069: u=gexp(gmul(pii2,z),prec);
1070: u1=gsub(u,gun);
1071: y=gdiv(gmul(gsqr(om2),elleisnum(om,2,0,prec)),pii2);
1072: y=gadd(ghalf,gdivgs(gmul(z,y),-12));
1073: y=gadd(y,ginv(u1));
1074: toadd=(long)ceil(9.065*gtodouble(gimag(z)));
1075: av1=avma; lim=stack_lim(av1,1); qn=q;
1076: do
1077: {
1078: p1=gadd(gdiv(u,gsub(gmul(qn,u),gun)),ginv(gsub(u,qn)));
1079: p1=gmul(qn,p1);
1080: y=gadd(y,p1);
1081: qn=gmul(q,qn);
1082: if (low_stack(lim, stack_lim(av1,1)))
1083: {
1084: GEN *gptr[2]; gptr[0]=&y; gptr[1]=&qn;
1085: if(DEBUGMEM>1) err(warnmem,"ellzeta");
1086: gerepilemany(av1,gptr,2);
1087: }
1088: }
1089: while (gexpo(qn) > - bit_accuracy(prec) - 5 - toadd);
1090:
1091: y=gmul(gdiv(pii2,om2),y);
1092: tetpil=avma;
1093: return gerepile(av,tetpil,gadd(y,et));
1094: }
1095:
1096: /* if flag=0, return ellsigma, otherwise return log(ellsigma) */
1097:
1098: static GEN
1099: ellsigmasum(GEN om, GEN z, long flag, long prec)
1100: {
1101: long av=avma,tetpil,lim,av1,si,toadd,n;
1102: GEN zinit,om1,om2,p1,pii2,tau,q,u,y,y1,qn,qn2,urn,urninv,v,ga;
1103: GEN uinv,x1,x2,et,etnew,uhalf,q8;
1104:
1105: if (!get_periods(om, &om1, &om2)) err(typeer,"ellsigmasum");
1106: p1=mppi(prec); setexpo(p1,2);
1107: pii2=cgetg(3,t_COMPLEX); pii2[1]=zero; pii2[2]=(long)p1;
1108: tau=gdiv(om1,om2); si=gsigne(gimag(tau));
1109: if (si==0)
1110: err(talker,"omega1 and omega2 are R-linearly dependent in ellsigma");
1111: if (si<0) { p1=om1; om1=om2; om2=p1; tau=ginv(tau); }
1112: v=getgamma(tau); tau=(GEN)v[2]; ga=(GEN)v[1];
1113: om2=gadd(gmul(gcoeff(ga,2,1),om1),gmul(gcoeff(ga,2,2),om2));
1114: om1=gmul(tau,om2); om=cgetg(3,t_VEC); om[1]=(long)om1; om[2]=(long)om2;
1115: z=gdiv(z,om2);
1116:
1117: x1=ground(gdiv(gimag(z),gimag(tau))); z=gsub(z,gmul(x1,tau));
1118: x2=ground(greal(z)); z=gsub(z,x2); zinit=gmul(z,om2);
1119: et=elleta(om,prec);
1120: etnew=gadd(gmul(x1,(GEN)et[1]),gmul(x2,(GEN)et[2]));
1121: etnew=gmul(etnew,gadd(gmul2n(gadd(gmul(x1,om1),gmul(x2,om2)),-1),zinit));
1122: if (mpodd(x1) || mpodd(x2)) etnew=gadd(etnew,gmul2n(pii2,-1));
1123: if (gexpo(z) < 5 - bit_accuracy(prec))
1124: {
1125: if (flag)
1126: {
1127: y=glog(zinit,prec);
1128: tetpil=avma;
1129: return gerepile(av,tetpil,gadd(etnew,y));
1130: }
1131: else
1132: {
1133: et=gexp(et,prec);
1134: tetpil=avma;
1135: return gerepile(av,tetpil,gmul(etnew,zinit));
1136: }
1137: }
1138:
1139: y1=gadd(etnew,gmul2n(gmul(gmul(z,zinit),(GEN)et[2]),-1));
1140:
1141: q8=gexp(gmul2n(gmul(pii2,tau),-3),prec);
1142: q=gpuigs(q8,8);
1143: uhalf=gexp(gmul(gmul2n(pii2,-1),z),prec);
1144: u=gneg_i(gsqr(uhalf)); uinv=ginv(u);
1145: y=gzero;
1146: toadd=(long)ceil(9.065*gtodouble(gabs(gimag(z),prec)));
1147: /* 9.065 = 2*Pi/log(2) */
1148: av1=avma; lim=stack_lim(av1,1); qn=q; qn2=gun;
1149: urn=uhalf; urninv=ginv(uhalf); n=0;
1150: do
1151: {
1152: y=gadd(y,gmul(qn2,gsub(urn,urninv)));
1153: qn2=gmul(qn,qn2);
1154: qn=gmul(q,qn);
1155: urn=gmul(urn,u); urninv=gmul(urninv,uinv);
1156: n++;
1157: if (low_stack(lim, stack_lim(av1,1)))
1158: {
1159: GEN *gptr[5]; gptr[0]=&y; gptr[1]=&qn; gptr[2]=&qn2; gptr[3]=&urn;
1160: gptr[4]=&urninv;
1161: if(DEBUGMEM>1) err(warnmem,"ellsigma");
1162: gerepilemany(av1,gptr,5);
1163: }
1164: }
1165: while (gexpo(qn2) + (n-1)*toadd > - bit_accuracy(prec) - 5);
1166:
1167: p1=gmul(q8,gmul(gdiv(gdiv((GEN)om[2],pii2),gpuigs(trueeta(tau,prec),3)),y));
1168: if (flag)
1169: {
1170: p1=glog(p1,prec); tetpil=avma;
1171: return gerepile(av,tetpil,gadd(y1,p1));
1172: }
1173: else
1174: {
1175: y=gexp(y1,prec); tetpil=avma;
1176: return gerepile(av,tetpil,gmul(p1,y));
1177: }
1178: }
1179:
1180: /* if flag=0, return ellsigma, otherwise return log(ellsigma) */
1181:
1182: static GEN
1183: ellsigmaprod(GEN om, GEN z, long flag, long prec)
1184: {
1185: long av=avma,tetpil,lim,av1,si,toadd;
1186: GEN zinit,om1,om2,p1,pii2,tau,q,u,y,y1,u1,qn,v,ga,negu,uinv,x1,x2,et,etnew,uhalf;
1187:
1188: if (!get_periods(om, &om1, &om2)) err(typeer,"ellsigmaprod");
1189: p1=mppi(prec); setexpo(p1,2);
1190: pii2=cgetg(3,t_COMPLEX); pii2[1]=zero; pii2[2]=(long)p1;
1191: tau=gdiv(om1,om2); si=gsigne(gimag(tau));
1192: if (si==0)
1193: err(talker,"omega1 and omega2 are R-linearly dependent in ellsigma");
1194: if (si<0) { p1=om1; om1=om2; om2=p1; tau=ginv(tau); }
1195: v=getgamma(tau); tau=(GEN)v[2]; ga=(GEN)v[1];
1196: om2=gadd(gmul(gcoeff(ga,2,1),om1),gmul(gcoeff(ga,2,2),om2));
1197: om1=gmul(tau,om2); om=cgetg(3,t_VEC); om[1]=(long)om1; om[2]=(long)om2;
1198: z=gdiv(z,om2);
1199:
1200: x1=ground(gdiv(gimag(z),gimag(tau))); z=gsub(z,gmul(x1,tau));
1201: x2=ground(greal(z)); z=gsub(z,x2); zinit=gmul(z,om2);
1202: et=elleta(om,prec);
1203: etnew=gadd(gmul(x1,(GEN)et[1]),gmul(x2,(GEN)et[2]));
1204: etnew=gmul(etnew,gadd(gmul2n(gadd(gmul(x1,om1),gmul(x2,om2)),-1),zinit));
1205: if (mpodd(x1) || mpodd(x2)) etnew=gadd(etnew,gmul2n(pii2,-1));
1206: if (gexpo(z) < 5 - bit_accuracy(prec))
1207: {
1208: if (flag)
1209: {
1210: y=glog(zinit,prec);
1211: tetpil=avma;
1212: return gerepile(av,tetpil,gadd(etnew,y));
1213: }
1214: else
1215: {
1216: et=gexp(et,prec);
1217: tetpil=avma;
1218: return gerepile(av,tetpil,gmul(etnew,zinit));
1219: }
1220: }
1221:
1222: y1=gadd(etnew,gmul2n(gmul(gmul(z,zinit),(GEN)et[2]),-1));
1223:
1224: q=gexp(gmul(pii2,tau),prec);
1225: uhalf=gexp(gmul(gmul2n(pii2,-1),z),prec); u=gsqr(uhalf);
1226: uinv=ginv(u);
1227: u1=gsub(uhalf,ginv(uhalf));
1228: y=gdiv(gmul(om2,u1),pii2);
1229: toadd=(long)ceil(9.065*gtodouble(gabs(gimag(z),prec)));
1230: /* 9.065 = 2*Pi/log(2) */
1231: av1=avma; lim=stack_lim(av1,1); qn=q;
1232: negu=stoi(-1);
1233: do
1234: {
1235: p1=gmul(gadd(gmul(qn,u),negu),gadd(gmul(qn,uinv),negu));
1236: p1=gdiv(p1,gsqr(gadd(qn,negu)));
1237: y=gmul(y,p1);
1238: qn=gmul(q,qn);
1239: if (low_stack(lim, stack_lim(av1,1)))
1240: {
1241: GEN *gptr[2]; gptr[0]=&y; gptr[1]=&qn;
1242: if(DEBUGMEM>1) err(warnmem,"ellsigma");
1243: gerepilemany(av1,gptr,2);
1244: }
1245: }
1246: while (gexpo(qn) > - bit_accuracy(prec) - 5 - toadd);
1247:
1248: if (flag)
1249: {
1250: p1=glog(y,prec);
1251: tetpil=avma;
1252: return gerepile(av,tetpil,gadd(y1,p1));
1253: }
1254: else
1255: {
1256: p1=gexp(y1,prec);
1257: tetpil=avma;
1258: return gerepile(av,tetpil,gmul(p1,y));
1259: }
1260: }
1261:
1262: GEN
1263: ellsigma(GEN om, GEN z, long flag, long prec)
1264: {
1265: if (flag>=2) return ellsigmaprod(om,z,flag&1,prec);
1266: else return ellsigmasum(om,z,flag,prec);
1267: }
1268:
1269: GEN
1270: pointell(GEN e, GEN z, long prec)
1271: {
1272: long av=avma,tetpil;
1273: GEN y,yp,v,p1;
1274:
1275: checkbell(e);
1276: p1=weipellnumall((GEN)e[16],(GEN)e[15],z,1,prec);
1277: if (lg(p1)==2) { avma=av; v=cgetg(2,t_VEC); v[1]=zero; return v; }
1278: y = gsub((GEN)p1[1], gdivgs((GEN)e[6],12));
1279: yp = gsub((GEN)p1[2], gmul2n(ellLHS0(e,y),-1));
1280: tetpil=avma; v=cgetg(3,t_VEC); v[1]=lcopy(y); v[2]=lcopy(yp);
1281: return gerepile(av,tetpil,v);
1282: }
1283:
1284: GEN
1285: weipell(GEN e, long prec)
1286: {
1287: long av1,tetpil,precres,i,k,l;
1288: GEN res,p1,s,t;
1289:
1290: checkell(e); precres = 2*prec+2;
1291: res=cgetg(precres,t_SER);
1292: res[1] = evalsigne(1) | evalvalp(-2) | evalvarn(0);
1293: if (!prec) { setsigne(res,0); return res; }
1294: for (i=3; i<precres; i+=2) res[i]=zero;
1295: switch(prec)
1296: {
1297: default: res[8]=ldivgs((GEN)e[11],6048);
1298: case 3: res[6]=ldivgs((GEN)e[10],240);
1299: case 2: res[4]=zero;
1300: case 1: res[2]=un;
1301: case 0: break;
1302: }
1303: for (k=4; k<prec; k++)
1304: {
1305: av1 = avma;
1306: s = k&1? gzero: gsqr((GEN)res[k+2]);
1307: t = gzero;
1308: for (l=2; l+l<k; l++)
1309: t = gadd(t, gmul((GEN)res[(l+1)<<1],(GEN)res[(k-l+1)<<1]));
1310: p1=gmulsg(3,gadd(s,gmul2n(t,1)));
1311: tetpil=avma;
1312: p1=gdivgs(p1,(k-3)*(2*k+1));
1313: res[(k+1)<<1] = lpile(av1,tetpil,p1);
1314: }
1315: return res;
1316: }
1317:
1318: GEN
1319: ellwp0(GEN om, GEN z, long flag, long prec, long PREC)
1320: {
1321: GEN v,om1,om2;
1322: long av = avma;
1323:
1324: if (z==NULL) return weipell(om,PREC);
1325: if (typ(z)==t_POL)
1326: {
1327: if (lgef(z) != 4 || !gcmp0((GEN)z[2]) || !gcmp1((GEN)z[3]))
1328: err(talker,"expecting a simple variable in ellwp");
1329: v = weipell(om,PREC); setvarn(v, varn(z));
1330: return v;
1331: }
1332: if (!get_periods(om, &om1, &om2)) err(typeer,"ellwp");
1333: switch(flag)
1334: {
1335: case 0: v=weipellnumall(om1,om2,z,0,prec);
1336: if (typ(v)==t_VEC && lg(v)==2) { avma=av; v=gpuigs(z,-2); }
1337: return v;
1338: case 1: v=weipellnumall(om1,om2,z,1,prec);
1339: if (typ(v)==t_VEC && lg(v)==2)
1340: {
1341: GEN p1 = gmul2n(gpuigs(z,3),1);
1342: long tetpil=avma;
1343: v=cgetg(3,t_VEC);
1344: v[1]=lpuigs(z,-2);
1345: v[2]=lneg(p1); return gerepile(av,tetpil,v);
1346: }
1347: return v;
1348: case 2: return pointell(om,z,prec);
1349: default: err(flagerr,"ellwp"); return NULL;
1350: }
1351: }
1352:
1353: /* compute a_2 using Jacobi sum */
1354: static GEN
1355: _a_2(GEN e)
1356: {
1357: long av = avma;
1358: GEN unmodp = gmodulss(1,8);
1359: ulong e6 = itos((GEN)gmul(unmodp,(GEN)e[6])[2]);
1360: ulong e8 = itos((GEN)gmul(unmodp,(GEN)e[8])[2]);
1361: ulong e72= itos((GEN)gmul(unmodp,gmul2n((GEN)e[7],1))[2]);
1362: long s = kross(e8, 2) + kross(e8 + e72 + e6 + 4, 2);
1363: avma = av; return stoi(-s);
1364: }
1365:
1366: /* a_p using Jacobi sums */
1367: static GEN
1368: apell2_intern(GEN e, ulong p)
1369: {
1370: if (p == 2) return _a_2(e);
1371: else
1372: {
1373: long av=avma,i;
1374: GEN unmodp = gmodulss(1,p);
1375: ulong e6 = itos((GEN)gmul(unmodp,(GEN)e[6])[2]);
1376: ulong e8 = itos((GEN)gmul(unmodp,(GEN)e[8])[2]);
1377: ulong e72= itos((GEN)gmul(unmodp,(GEN)e[7])[2]) << 1;
1378: long s = kross(e8, p);
1379:
1380: if (p < 757)
1381: for (i=1; i<p; i++)
1382: s += kross(e8 + i*(e72 + i*(e6 + (i<<2))), p);
1383: else
1384: for (i=1; i<p; i++)
1385: s += kross(e8 + mulssmod(i, e72 + mulssmod(i, e6 + (i<<2), p), p), p);
1386: avma=av; return stoi(-s);
1387: }
1388: }
1389:
1390: GEN
1391: apell2(GEN e, GEN pp)
1392: {
1393: checkell(e); if (typ(pp)!=t_INT) err(elliper1);
1394: if (expi(pp) > 29)
1395: err(talker,"prime too large in jacobi apell2, use apell instead");
1396:
1397: return apell2_intern(e, (ulong)pp[2]);
1398: }
1399:
1400: GEN ellap0(GEN e, GEN p, long flag)
1401: {
1402: return flag? apell2(e,p): apell(e,p);
1403: }
1404:
1405: /* invert all elements of x mod p using Montgomery's trick */
1406: GEN
1407: multi_invmod(GEN x, GEN p)
1408: {
1409: long i, lx = lg(x);
1410: GEN u,y = cgetg(lx, t_VEC);
1411:
1412: y[1] = x[1];
1413: for (i=2; i<lx; i++)
1414: y[i] = lresii(mulii((GEN)y[i-1], (GEN)x[i]), p);
1415:
1416: u = mpinvmod((GEN)y[--i], p);
1417: for ( ; i > 1; i--)
1418: {
1419: y[i] = lresii(mulii(u, (GEN)y[i-1]), p);
1420: u = resii(mulii(u, (GEN)x[i]), p); /* u = 1 / (x[1] ... x[i-1]) */
1421: }
1422: y[1] = (long)u; return y;
1423: }
1424:
1425: static GEN
1426: addsell(GEN e, GEN z1, GEN z2, GEN p)
1427: {
1428: GEN p1,p2,x,x1,x2,y,y1,y2;
1429: long av = avma;
1430:
1431: if (!z1) return z2;
1432: if (!z2) return z1;
1433:
1434: x1 = (GEN)z1[1]; y1 = (GEN)z1[2];
1435: x2 = (GEN)z2[1]; y2 = (GEN)z2[2];
1436: p2 = subii(x2, x1);
1437: if (p2 == gzero)
1438: {
1439: if (!signe(y1) || !egalii(y1,y2)) return NULL;
1440: p2 = shifti(y1,1);
1441: p1 = addii(e, mulii(x1,mulsi(3,x1)));
1442: p1 = resii(p1, p);
1443: }
1444: else p1 = subii(y2,y1);
1445: p1 = mulii(p1, mpinvmod(p2, p));
1446: p1 = resii(p1, p);
1447: x = subii(sqri(p1), addii(x1,x2)); x = modii(x,p);
1448: y = negi(addii(y1, mulii(p1,subii(x,x1))));
1449: avma = av; p1 = cgetg(3,t_VEC);
1450: p1[1] = licopy(x);
1451: p1[2] = lmodii(y, p); return p1;
1452: }
1453:
1454: /* z1 <-- z1 + z2 */
1455: static void
1456: addsell_part2(GEN e, GEN z1, GEN z2, GEN p, GEN p2inv)
1457: {
1458: GEN p1,x,x1,x2,y,y1,y2;
1459:
1460: x1 = (GEN)z1[1]; y1 = (GEN)z1[2];
1461: x2 = (GEN)z2[1]; y2 = (GEN)z2[2];
1462: if (x1 == x2)
1463: {
1464: p1 = addii(e, mulii(x1,mulsi(3,x1)));
1465: p1 = resii(p1, p);
1466: }
1467: else p1 = subii(y2,y1);
1468:
1469: p1 = mulii(p1, p2inv);
1470: p1 = resii(p1, p);
1471: x = subii(sqri(p1), addii(x1,x2)); x = modii(x,p);
1472: y = negi(addii(y1, mulii(p1,subii(x,x1)))); y = modii(y,p);
1473: affii(x, x1);
1474: affii(y, y1);
1475: }
1476:
1477: static GEN
1478: powsell(GEN e, GEN z, GEN n, GEN p)
1479: {
1480: GEN y;
1481: long s=signe(n),i,j;
1482: ulong m;
1483:
1484: if (!s || !z) return NULL;
1485: if (s < 0)
1486: {
1487: n = negi(n); y = cgetg(3,t_VEC);
1488: y[2] = lnegi((GEN)z[2]);
1489: y[1] = z[1]; z = y;
1490: }
1491: if (is_pm1(n)) return z;
1492:
1493: y = NULL;
1494: for (i=lgefint(n)-1; i>2; i--)
1495: for (m=n[i],j=0; j<BITS_IN_LONG; j++,m>>=1)
1496: {
1497: if (m&1) y = addsell(e,y,z,p);
1498: z = addsell(e,z,z,p);
1499: }
1500: for (m=n[2]; m>1; m>>=1)
1501: {
1502: if (m&1) y = addsell(e,y,z,p);
1503: z = addsell(e,z,z,p);
1504: }
1505: return addsell(e,y,z,p);
1506: }
1507:
1508: /* make sure *x has lgefint >= k */
1509: static void
1510: _fix(GEN x, long k)
1511: {
1512: GEN y = (GEN)*x;
1513: if (lgefint(y) < k) { GEN p1 = cgeti(k); affii(y,p1); *x = (long)p1; }
1514: }
1515:
1516: /* low word of integer x */
1517: #define _low(x) ((x)[lgefint(x)-1])
1518:
1519: /* compute a_p using Shanks/Mestre + Montgomery's trick. Assume p > 20, say */
1520: GEN
1521: apell1(GEN e, GEN p)
1522: {
1523: long *tx, *ty, *ti, av = avma, av2,pfinal,i,j,j2,s,flc,flcc,x,nb;
1524: GEN p1,p2,p3,h,mfh,f,fh,fg,pordmin,u,v,p1p,p2p,acon,bcon,c4,c6,cp4,pts;
1525:
1526: if (DEBUGLEVEL) timer2();
1527: p1 = gmodulsg(1,p);
1528: c4 = gdivgs(gmul(p1,(GEN)e[10]), -48);
1529: c6 = gdivgs(gmul(p1,(GEN)e[11]), -864);
1530: pordmin = gceil(gmul2n(gsqrt(p,DEFAULTPREC),2));
1531: p1p = addsi(1,p); p2p = shifti(p1p,1);
1532: x=0; flcc=0; flc = kronecker((GEN)c6[2],p);
1533: u=c6; acon=gzero; bcon=gun; h=p1p;
1534: for(;;)
1535: {
1536: while (flc==flcc || !flc)
1537: {
1538: x++;
1539: u = gadd(c6, gmulsg(x, gaddgs(c4,x*x)));
1540: flc = kronecker((GEN)u[2],p);
1541: }
1542: flcc = flc;
1543: f = cgetg(3,t_VEC);
1544: f[1] = (long)lift_intern(gmulsg(x,u));
1545: f[2] = (long)lift_intern(gsqr(u));
1546: cp4 = lift_intern(gmul(c4, (GEN)f[2]));
1547: fh = powsell(cp4,f,h,p);
1548: s = itos(gceil(gsqrt(gdiv(pordmin,bcon),DEFAULTPREC))) >> 1;
1549: nb = min(128, s >> 1);
1550: if (bcon == gun)
1551: { /* first time: initialize */
1552: tx = newbloc(s+1); *tx = evaltyp(t_VECSMALL) | evallg(s+1);
1553: ty = newbloc(s+1);
1554: ti = newbloc(s+1);
1555: }
1556: else f = powsell(cp4,f,bcon,p);
1557: if (!fh) goto FOUND;
1558:
1559: p1 = gcopy(fh);
1560: pts = new_chunk(nb+1);
1561: j = lgefint(p);
1562: for (i=1; i<=nb; i++)
1563: { /* baby steps */
1564: pts[i] = (long)p1;
1565: _fix(p1+1, j); tx[i] = _low((GEN)p1[1]);
1566: _fix(p1+2, j); ty[i] = _low((GEN)p1[2]);
1567: p1 = addsell(cp4,p1,f,p); /* f^h * F^nb */
1568: if (!p1) { h = addii(h, mulsi(i,bcon)); goto FOUND; }
1569: }
1570: mfh = dummycopy(fh);
1571: mfh[2] = lnegi((GEN)mfh[2]);
1572: fg = addsell(cp4,p1,mfh,p); /* F^nb */
1573: if (!fg) { h = addii(h, mulsi(nb,bcon)); goto FOUND; }
1574: u = cgetg(nb+1, t_VEC);
1575: av2 = avma; /* more baby steps, nb points at a time */
1576: while (i <= s)
1577: {
1578: long maxj;
1579: for (j=1; j<=nb; j++)
1580: {
1581: p1 = (GEN)pts[j];
1582: u[j] = lsubii((GEN)fg[1], (GEN)p1[1]);
1583: if (u[j] == zero)
1584: {
1585: if (!signe(p1[2]) || !egalii((GEN)p1[2],(GEN)fg[2]))
1586: { h = addii(h, mulsi(i+j-1,bcon)); goto FOUND; }
1587: /* doubling never occurs */
1588: err(bugparier, "apell1: doubling?");
1589: }
1590: }
1591: v = multi_invmod(u, p);
1592: maxj = (i-1 + nb <= s)? nb: s % nb;
1593: for (j=1; j<=maxj; j++,i++)
1594: {
1595: p1 = (GEN)pts[j];
1596: addsell_part2(cp4,p1,fg,p, (GEN)v[j]);
1597: tx[i] = _low((GEN)p1[1]);
1598: ty[i] = _low((GEN)p1[2]);
1599: }
1600: avma = av2;
1601: }
1602: p1 = addsell(cp4,(GEN)pts[j-1],mfh,p); /* = f^(s-1) */
1603: if (DEBUGLEVEL) msgtimer("[apell1] baby steps, s = %ld",s);
1604:
1605: /* giant steps: fg = f^s */
1606: fg = addsell(cp4,p1,f,p);
1607: if (!fg) { h = addii(h, mulsi(s,bcon)); goto FOUND; }
1608: pfinal = _low(p); av2 = avma;
1609:
1610: p1 = gen_sort(tx, cmp_IND | cmp_C, NULL);
1611: for (i=1; i<=s; i++) ti[i] = tx[p1[i]];
1612: for (i=1; i<=s; i++) { tx[i] = ti[i]; ti[i] = ty[p1[i]]; }
1613: for (i=1; i<=s; i++) { ty[i] = ti[i]; ti[i] = p1[i]; }
1614: if (DEBUGLEVEL) msgtimer("[apell1] sorting");
1615: avma = av2;
1616:
1617: gaffect(fg, (GEN)pts[1]);
1618: for (j=2; j<=nb; j++) /* pts = first nb multiples of fg */
1619: gaffect(addsell(cp4,(GEN)pts[j-1],fg,p), (GEN)pts[j]);
1620: /* replace fg by fg^nb since we do nb points at a time */
1621: avma = av2;
1622: fg = gcopy((GEN)pts[nb]);
1623: av2 = avma;
1624:
1625: for (i=1,j=1; ; i++)
1626: {
1627: GEN ftest = (GEN)pts[j];
1628: ulong m, l = 1, r = s+1;
1629: long k, k2;
1630:
1631: avma = av2;
1632: k = _low((GEN)ftest[1]);
1633: while (l<r)
1634: {
1635: m = (l+r) >> 1;
1636: if (tx[m] < k) l = m+1; else r = m;
1637: }
1638: if (r <= s && tx[r] == k)
1639: {
1640: while (tx[r] == k && r) r--;
1641: k2 = _low((GEN)ftest[2]);
1642: for (r++; tx[r] == k && r <= s; r++)
1643: if (ty[r] == k2 || ty[r] == pfinal - k2)
1644: { /* [h+j2] f == ± ftest (= [i.s] f)? */
1645: if (DEBUGLEVEL) msgtimer("[apell1] giant steps, i = %ld",i);
1646: j2 = ti[r] - 1;
1647: p1 = addsell(cp4, powsell(cp4,f,stoi(j2),p),fh,p);
1648: if (egalii((GEN)p1[1], (GEN)ftest[1]))
1649: {
1650: h = addii(h, mulsi(j2,bcon));
1651: if (egalii((GEN)p1[2], (GEN)ftest[2])) i = -i;
1652: h = addii(h, mulsi(s, mulsi(i, bcon)));
1653: goto FOUND;
1654: }
1655: }
1656: }
1657: if (++j > nb)
1658: { /* compute next nb points */
1659: long save;
1660: for (j=1; j<=nb; j++)
1661: {
1662: p1 = (GEN)pts[j];
1663: u[j] = lsubii((GEN)fg[1], (GEN)p1[1]);
1664: if (u[j] == zero) /* occurs once: i = j = nb, p1 == fg */
1665: {
1666: u[j] = lshifti((GEN)p1[2],1);
1667: save = fg[1]; fg[1] = p1[1];
1668: }
1669: }
1670: v = multi_invmod(u, p);
1671: for (j=1; j<=nb; j++)
1672: addsell_part2(cp4, (GEN)pts[j],fg,p, (GEN)v[j]);
1673: if (i == nb) { fg[1] = save; }
1674: j = 1;
1675: }
1676: }
1677:
1678: FOUND: /* success, found a point of exponent h */
1679: p2 = decomp(h); p1=(GEN)p2[1]; p2=(GEN)p2[2];
1680: for (i=1; i<lg(p1); i++)
1681: for (j=itos((GEN)p2[i]); j; j--)
1682: {
1683: p3 = divii(h,(GEN)p1[i]);
1684: if (powsell(cp4,f,p3,p)) break;
1685: h = p3;
1686: }
1687: /* now h is the exact order */
1688: if (bcon == gun) bcon = h;
1689: else
1690: {
1691: p1 = chinois(gmodulcp(acon,bcon), gmodulsg(0,h));
1692: acon = (GEN)p1[2];
1693: bcon = (GEN)p1[1];
1694: }
1695:
1696: i = (cmpii(bcon, pordmin) < 0);
1697: if (i) acon = centermod(subii(p2p,acon), bcon);
1698: p1 = ground(gdiv(gsub(p1p,acon),bcon));
1699: h = addii(acon, mulii(p1,bcon));
1700: if (!i) break;
1701: }
1702: gunclone(tx);
1703: gunclone(ty);
1704: gunclone(ti);
1705: p1 = (flc==1)? subii(p1p,h): subii(h,p1p);
1706: return gerepileupto(av,p1);
1707: }
1708:
1709: typedef struct
1710: {
1711: int isnull;
1712: long x,y;
1713: } sellpt;
1714:
1715: /* set p1 <-- p1 + p2, safe with p1 = p2 */
1716: static void
1717: addsell1(long e, long p, sellpt *p1, sellpt *p2)
1718: {
1719: long num, den, lambda;
1720:
1721: if (p1->isnull) { *p1 = *p2; return; }
1722: if (p2->isnull) return;
1723: if (p1->x == p2->x)
1724: {
1725: if (! p1->y || p1->y != p2->y) { p1->isnull = 1; return; }
1726: num = addssmod(e, mulssmod(3, mulssmod(p1->x, p1->x, p), p), p);
1727: den = addssmod(p1->y, p1->y, p);
1728: }
1729: else
1730: {
1731: num = subssmod(p1->y, p2->y, p);
1732: den = subssmod(p1->x, p2->x, p);
1733: }
1734: lambda = divssmod(num, den, p);
1735: num = subssmod(mulssmod(lambda, lambda, p), addssmod(p1->x, p2->x, p), p);
1736: p1->y = subssmod(mulssmod(lambda, subssmod(p2->x, num, p), p), p2->y, p);
1737: p1->x = num; /* necessary in case p1 = p2: we need p2->x above */
1738: }
1739:
1740: static void
1741: powssell1(long e, long p, long n, sellpt *p1, sellpt *p2)
1742: {
1743: sellpt p3 = *p1;
1744:
1745: if (n < 0) { n = -n; if (p3.y) p3.y = p - p3.y; }
1746: p2->isnull = 1;
1747: for(;;)
1748: {
1749: if (n&1) addsell1(e, p, p2, &p3);
1750: n>>=1; if (!n) return;
1751: addsell1(e, p, &p3, &p3);
1752: }
1753: }
1754:
1755: typedef struct
1756: {
1757: long x,y,i;
1758: } multiple;
1759:
1760: static int
1761: compare_multiples(multiple *a, multiple *b)
1762: {
1763: return a->x - b->x;
1764: }
1765:
1766: /* assume e has good reduction at p. Should use Montgomery. */
1767: static GEN
1768: apell0(GEN e, long p)
1769: {
1770: GEN p1,p2;
1771: sellpt f,fh,fg,ftest,f2;
1772: long pordmin,u,p1p,p2p,acon,bcon,c4,c6,cp4;
1773: long av,i,j,s,flc,flcc,x,q,h,p3,l,r,m;
1774: multiple *table;
1775:
1776: if (p < 99) return apell2_intern(e,(ulong)p);
1777:
1778: av = avma; p1 = gmodulss(1,p);
1779: c4 = itos((GEN)gdivgs(gmul(p1,(GEN)e[10]), -48)[2]);
1780: c6 = itos((GEN)gdivgs(gmul(p1,(GEN)e[11]), -864)[2]);
1781: pordmin = (long)(1 + 4*sqrt((float)p));
1782: p1p = p+1; p2p = p1p << 1;
1783: x=0; flcc=0; flc = kross(c6, p);
1784: u=c6; acon=0; bcon=1; h=p1p;
1785: for(;;)
1786: {
1787: while (flc==flcc || !flc)
1788: {
1789: x++;
1790: u = addssmod(c6, mulssmod(x, c4+mulssmod(x,x,p), p), p);
1791: flc = kross(u,p);
1792: }
1793: flcc = flc;
1794: f.isnull = 0;
1795: f.x = mulssmod(x, u, p);
1796: f.y = mulssmod(u, u, p);
1797: cp4 = mulssmod(c4, f.y, p);
1798: powssell1(cp4, p, h, &f, &fh);
1799: s = (long) (sqrt(((float)pordmin)/bcon) / 2);
1800: if (!s) s=1;
1801: if (bcon==1)
1802: {
1803: table = (multiple *) gpmalloc((s+1)*sizeof(multiple));
1804: f2 = f;
1805: }
1806: else powssell1(cp4, p, bcon, &f, &f2);
1807: for (i=0; i < s; i++)
1808: {
1809: if (fh.isnull) { h += bcon*i; goto FOUND; }
1810: table[i].x = fh.x;
1811: table[i].y = fh.y;
1812: table[i].i = i;
1813: addsell1(cp4, p, &fh, &f2);
1814: }
1815: qsort(table,s,sizeof(multiple),(QSCOMP)compare_multiples);
1816: powssell1(cp4, p, s, &f2, &fg); ftest = fg;
1817: for (i=1; ; i++)
1818: {
1819: if (ftest.isnull) err(bugparier,"apell (f^(i*s) = 1)");
1820: l=0; r=s;
1821: while (l<r)
1822: {
1823: m = (l+r) >> 1;
1824: if (table[m].x < ftest.x) l=m+1; else r=m;
1825: }
1826: if (r < s && table[r].x == ftest.x) break;
1827: addsell1(cp4, p, &ftest, &fg);
1828: }
1829: h += table[r].i * bcon;
1830: if (table[r].y == ftest.y) i = -i;
1831: h += s * i * bcon;
1832:
1833: FOUND:
1834: p2=decomp(stoi(h)); p1=(GEN)p2[1]; p2=(GEN)p2[2];
1835: for (i=1; i < lg(p1); i++)
1836: for (j = mael(p2,i,2); j; j--)
1837: {
1838: p3 = h / mael(p1,i,2);
1839: powssell1(cp4, p, p3, &f, &fh);
1840: if (!fh.isnull) break;
1841: h = p3;
1842: }
1843: if (bcon == 1) bcon = h;
1844: else
1845: {
1846: p1 = chinois(gmodulss(acon,bcon), gmodulss(0,h));
1847: acon = itos((GEN)p1[2]);
1848: if (is_bigint(p1[1])) { h = acon; break; }
1849: bcon = itos((GEN)p1[1]);
1850: }
1851:
1852: i = (bcon < pordmin);
1853: if (i)
1854: {
1855: acon = (p2p - acon) % bcon;
1856: if ((acon << 1) > bcon) acon -= bcon;
1857: }
1858: q = ((ulong)(p2p + bcon - (acon << 1))) / (bcon << 1);
1859: h = acon + q*bcon;
1860: avma = av; if (!i) break;
1861: }
1862: free(table); return stoi((flc==1)? p1p-h: h-p1p);
1863: }
1864:
1865: GEN
1866: apell(GEN e, GEN p)
1867: {
1868: checkell(e);
1869: if (typ(p)!=t_INT || signe(p)<0) err(talker,"not a prime in apell");
1870: if (gdivise((GEN)e[12],p)) /* e[12] may be an intmod */
1871: {
1872: long av = avma,s;
1873: GEN c6 = gmul((GEN)e[11],gmodulsg(1,p));
1874: s = kronecker((GEN)c6[2],p); avma=av;
1875: switch(mod4(p))
1876: {
1877: case 0:
1878: case 3: s = -s;
1879: }
1880: return stoi(s);
1881: }
1882: if (cmpis(p, 0x3fffffff) > 0) return apell1(e, p);
1883: return apell0(e, itos(p));
1884: }
1885:
1886: /* TEMPC is the largest prime whose square is less than HIGHBIT */
1887: #ifndef LONG_IS_64BIT
1888: # define TEMPC 46337
1889: # define TEMPMAX 16777215UL
1890: #else
1891: # define TEMPC 3037000493
1892: # define TEMPMAX 4294967295UL
1893: #endif
1894:
1895: GEN
1896: anell(GEN e, long n)
1897: {
1898: long tab[4]={0,1,1,-1}; /* p prime; (-1/p) = tab[p&3]. tab[0] is not used */
1899: long p, i, m, av, tetpil;
1900: GEN p1,p2,an;
1901:
1902: checkell(e);
1903: if (n <= 0) return cgetg(1,t_VEC);
1904: if (n>TEMPMAX) err(impl,"anell for n>=2^24 (or 2^32 for 64 bit machines)");
1905: an = cgetg(n+1,t_VEC); an[1] = un;
1906: for (i=2; i <= n; i++) an[i] = 0;
1907: for (p=2; p<=n; p++)
1908: if (!an[p])
1909: {
1910: if (!smodis((GEN)e[12],p)) /* mauvaise reduction, p | e[12] */
1911: switch (tab[p&3] * krogs((GEN)e[11],p)) /* renvoie (-c6/p) */
1912: {
1913: case -1: /* non deployee */
1914: for (m=p; m<=n; m+=p)
1915: if (an[m/p]) an[m]=lneg((GEN)an[m/p]);
1916: continue;
1917: case 0: /* additive */
1918: for (m=p; m<=n; m+=p) an[m]=zero;
1919: continue;
1920: case 1: /* deployee */
1921: for (m=p; m<=n; m+=p)
1922: if (an[m/p]) an[m]=lcopy((GEN)an[m/p]);
1923: }
1924: else /* bonne reduction */
1925: {
1926: GEN ap = apell0(e,p);
1927:
1928: if (p < TEMPC)
1929: {
1930: ulong pk, oldpk = 1;
1931: for (pk=p; pk <= n; oldpk=pk, pk *= p)
1932: {
1933: if (pk == p) an[pk] = (long) ap;
1934: else
1935: {
1936: av = avma;
1937: p1 = mulii(ap, (GEN)an[oldpk]);
1938: p2 = mulsi(p, (GEN)an[oldpk/p]);
1939: tetpil = avma;
1940: an[pk] = lpile(av,tetpil,subii(p1,p2));
1941: }
1942: for (m = n/pk; m > 1; m--)
1943: if (an[m] && m%p) an[m*pk] = lmulii((GEN)an[m], (GEN)an[pk]);
1944: }
1945: }
1946: else
1947: {
1948: an[p] = (long) ap;
1949: for (m = n/p; m > 1; m--)
1950: if (an[m] && m%p) an[m*p] = lmulii((GEN)an[m], (GEN)an[p]);
1951: }
1952: }
1953: }
1954: return an;
1955: }
1956:
1957: GEN
1958: akell(GEN e, GEN n)
1959: {
1960: long i,j,ex,av=avma;
1961: GEN p1,p2,ap,u,v,w,y,pl;
1962:
1963: checkell(e);
1964: if (typ(n)!=t_INT) err(talker,"not an integer type in akell");
1965: if (signe(n)<= 0) return gzero;
1966: y=gun; if (gcmp1(n)) return y;
1967: p2=auxdecomp(n,1); p1=(GEN)p2[1]; p2=(GEN)p2[2];
1968: for (i=1; i<lg(p1); i++)
1969: {
1970: pl=(GEN)p1[i];
1971: if (divise((GEN)e[12], pl)) /* mauvaise reduction */
1972: {
1973: j = (((mod4(pl)+1)&2)-1)*kronecker((GEN)e[11],pl);
1974: if (j<0 && mpodd((GEN)p2[i])) y = negi(y);
1975: if (!j) { avma=av; return gzero; }
1976: }
1977: else /* bonne reduction */
1978: {
1979: ap=apell(e,pl); ex=itos((GEN)p2[i]);
1980: u=ap; v=gun;
1981: for (j=2; j<=ex; j++)
1982: {
1983: w = subii(mulii(ap,u), mulii(pl,v));
1984: v=u; u=w;
1985: }
1986: y = mulii(u,y);
1987: }
1988: }
1989: return gerepileupto(av,y);
1990: }
1991:
1992: GEN
1993: hell(GEN e, GEN a, long prec)
1994: {
1995: long av=avma,tetpil,n;
1996: GEN p1,p2,y,z,q,pi2surw,pi2isurw,qn,ps;
1997:
1998: checkbell(e);
1999: pi2surw=gdiv(gmul2n(mppi(prec),1),(GEN)e[15]);
2000: pi2isurw=cgetg(3,t_COMPLEX); pi2isurw[1]=zero; pi2isurw[2]=(long)pi2surw;
2001: z=gmul(greal(zell(e,a,prec)),pi2surw);
2002: q=greal(gexp(gmul((GEN)e[16],pi2isurw),prec));
2003: y=gsin(z,prec); n=0; qn=gun; ps=gneg_i(q);
2004: do
2005: {
2006: n++; p1=gsin(gmulsg(2*n+1,z),prec); qn=gmul(qn,ps);
2007: ps=gmul(ps,q); p1=gmul(p1,qn); y=gadd(y,p1);
2008: }
2009: while (gexpo(qn) >= - bit_accuracy(prec));
2010: p1=gmul(gsqr(gdiv(gmul2n(y,1), d_ellLHS(e,a))),pi2surw);
2011: p2=gsqr(gsqr(gdiv(p1,gsqr(gsqr(denom((GEN)a[1]))))));
2012: p1=gdiv(gmul(p2,q),(GEN)e[12]);
2013: p1=gmul2n(glog(gabs(p1,prec),prec),-5);
2014: tetpil=avma; return gerepile(av,tetpil,gneg(p1));
2015: }
2016:
2017: static GEN
2018: hells(GEN e, GEN x, long prec)
2019: {
2020: GEN w,z,t,mu,e72,e82;
2021: long n,lim;
2022:
2023: t = gdiv(realun(prec),(GEN)x[1]);
2024: mu = gmul2n(glog(numer((GEN)x[1]),prec),-1);
2025: e72 = gmul2n((GEN)e[7],1);
2026: e82 = gmul2n((GEN)e[8],1);
2027: lim = 6 + (bit_accuracy(prec) >> 1);
2028: for (n=0; n<lim; n++)
2029: {
2030: w = gmul(t,gaddsg(4,gmul(t,gadd((GEN)e[6],gmul(t,gadd(e72,gmul(t,(GEN)e[8])))))));
2031: z = gsub(gun,gmul(gsqr(t),gadd((GEN)e[7],gmul(t,gadd(e82,gmul(t,(GEN)e[9]))))));
2032: mu = gadd(mu,gmul2n(glog(z,prec), -((n<<1)+3)));
2033: t = gdiv(w,z);
2034: }
2035: return mu;
2036: }
2037:
2038: GEN
2039: hell2(GEN e, GEN x, long prec)
2040: {
2041: GEN ep,e3,ro,p1,p2,mu,d,xp;
2042: long av=avma,tetpil,lx,lc,i,j,tx;
2043:
2044: if (!oncurve(e,x)) err(heller1);
2045: d=(GEN)e[12]; ro=(GEN)e[14]; e3=(gsigne(d) < 0)?(GEN)ro[1]:(GEN)ro[3];
2046: p1=cgetg(5,t_VEC); p1[1]=un; p1[2]=laddgs(gfloor(e3),-1);
2047: p1[3]=p1[4]=zero; ep=coordch(e,p1); xp=pointch(x,p1);
2048: tx=typ(x[1]); lx=lg(x);
2049: if (!is_matvec_t(tx))
2050: {
2051: if (lx<3) { avma=av; return gzero; }
2052: tetpil=avma; return gerepile(av,tetpil,hells(ep,xp,prec));
2053: }
2054: tx=typ(x);
2055: tetpil=avma; mu=cgetg(lx,tx);
2056: if (tx != t_MAT)
2057: for (i=1; i<lx; i++) mu[i]=(long)hells(ep,(GEN)xp[i],prec);
2058: else
2059: {
2060: lc=lg(x[1]);
2061: for (i=1; i<lx; i++)
2062: {
2063: p1=cgetg(lc,t_COL); mu[i]=(long)p1; p2=(GEN)xp[i];
2064: for (j=1; j<lc; j++) p1[j]=(long)hells(ep,(GEN)p2[j],prec);
2065: }
2066: }
2067: return gerepile(av,tetpil,mu);
2068: }
2069:
2070: GEN
2071: hell0(GEN e, GEN z, long prec)
2072: {
2073: GEN a,b,s,x,u,v,u1,p1,p2;
2074: long n,i, ex = 5-bit_accuracy(prec);
2075:
2076: /* cf. zell mais ne marche pas. Comment corriger? K.B. */
2077: x = new_coords(e,(GEN)z[1],&a,&b,prec);
2078:
2079: u = gmul2n(gadd(a,b), -1);
2080: v = gsqrt(gmul(a,b), prec); s = gun;
2081: for(n=0; ; n++)
2082: {
2083: p1 = gmul2n(gsub(x, gsqr(v)), -1);
2084: p2 = gsqr(u);
2085: x = gadd(p1, gsqrt(gadd(gsqr(p1), gmul(x, p2)), prec));
2086: p2 = gadd(x, p2);
2087: for (i=1; i<=n; i++) p2 = gsqr(p2);
2088: s = gmul(s, p2);
2089: u1 = gmul2n(gadd(u,v), -1);
2090: if (gexpo(gsub(u,u1)) < ex) break;
2091:
2092: v = gsqrt(gmul(u,v), prec);
2093: u = u1;
2094: }
2095: return gmul2n(glog(gdiv(gsqr(p2), s), prec) ,-1);
2096: }
2097:
2098: /* On suppose que `e' est a coeffs entiers donnee par un modele minimal */
2099: static GEN
2100: ghell0(GEN e, GEN a, long flag, long prec)
2101: {
2102: long av=avma,lx,i,n,n2,grandn,tx=typ(a);
2103: GEN p,p1,p2,x,y,z,phi2,psi2,psi3,logdep;
2104:
2105: checkbell(e); if (!is_matvec_t(tx)) err(elliper1);
2106: lx = lg(a); if (lx==1) return cgetg(1,tx);
2107: tx=typ(a[1]);
2108: if (is_matvec_t(tx))
2109: {
2110: z=cgetg(lx,tx);
2111: for (i=1; i<lx; i++) z[i]=(long)ghell0(e,(GEN)a[i],flag,prec);
2112: return z;
2113: }
2114: if (lg(a)<3) return gzero;
2115: if (!oncurve(e,a)) err(heller1);
2116:
2117: psi2=numer(d_ellLHS(e,a));
2118: if (!signe(psi2)) { avma=av; return gzero; }
2119:
2120: x=(GEN)a[1]; y=(GEN)a[2];
2121: p2=gadd(gmulsg(3,(GEN)e[7]),gmul(x,gadd((GEN)e[6],gmulsg(3,x))));
2122: psi3=numer(gadd((GEN)e[9],gmul(x,gadd(gmulsg(3,(GEN)e[8]),gmul(x,p2)))));
2123: if (!signe(psi3)) { avma=av; return gzero; }
2124:
2125: p1 = gmul(x,gadd(shifti((GEN)e[2],1),gmulsg(3,x)));
2126: phi2=numer(gsub(gadd((GEN)e[4],p1), gmul((GEN)e[1],y)));
2127: p1=(GEN)factor(mppgcd(psi2,phi2))[1]; lx=lg(p1);
2128: switch(flag)
2129: {
2130: case 0: z = hell2(e,a,prec); break; /* Tate 4^n */
2131: case 1: z = hell(e,a,prec); break; /* Silverman's trick */
2132: case 2: z = hell0(e,a,prec); break; /* Mestre's trick */
2133:
2134: }
2135: for (i=1; i<lx; i++)
2136: {
2137: p=(GEN)p1[i];
2138: if (signe(resii((GEN)e[10],p)))
2139: {
2140: grandn=ggval((GEN)e[12],p);
2141: if (grandn)
2142: {
2143: n2=ggval(psi2,p); n=n2<<1;
2144: logdep=gneg_i(glog(p,prec));
2145: if (n>grandn) n=grandn;
2146: p2=divrs(mulsr(n*(grandn+grandn-n),logdep),grandn<<3);
2147: z=gadd(z,p2);
2148: }
2149: }
2150: else
2151: {
2152: n2=ggval(psi2,p);
2153: logdep=gneg_i(glog(p,prec));
2154: n=ggval(psi3,p);
2155: if (n>=3*n2) p2=gdivgs(mulsr(n2,logdep),3);
2156: else p2=gmul2n(mulsr(n,logdep),-3);
2157: z=gadd(z,p2);
2158: }
2159: }
2160: return gerepileupto(av,z);
2161: }
2162:
2163: GEN
2164: ellheight0(GEN e, GEN a, long flag, long prec)
2165: {
2166: switch(flag)
2167: {
2168: case 0: return ghell(e,a,prec);
2169: case 1: return ghell2(e,a,prec);
2170: case 2: return ghell0(e,a,2,prec);
2171: }
2172: err(flagerr,"ellheight");
2173: return NULL; /* not reached */
2174: }
2175:
2176: GEN
2177: ghell2(GEN e, GEN a, long prec)
2178: {
2179: return ghell0(e,a,0,prec);
2180: }
2181:
2182: GEN
2183: ghell(GEN e, GEN a, long prec)
2184: {
2185: return ghell0(e,a,1,prec);
2186: }
2187:
2188: static long ellrootno_all(GEN e, GEN p, GEN* ptcond);
2189:
2190: GEN
2191: lseriesell(GEN e, GEN s, GEN A, long prec)
2192: {
2193: long av=avma,av1,tetpil,lim,l,n,eps,flun;
2194: GEN z,p1,p2,cg,cg1,v,cga,cgb,s2,ns,gs,N;
2195:
2196: if (!A) A = gun;
2197: else
2198: {
2199: if (gsigne(A)<=0)
2200: err(talker,"cut-off point must be positive in lseriesell");
2201: if (gcmpgs(A,1) < 0) A = ginv(A);
2202: }
2203: flun = gcmp1(A) && gcmp1(s);
2204: eps = ellrootno_all(e,gun,&N);
2205: if (flun && eps<0) { z=cgetr(prec); affsr(0,z); return z; }
2206: cg1=mppi(prec); setexpo(cg1,2); cg=divrr(cg1,gsqrt(N,prec));
2207: cga=gmul(cg,A); cgb=gdiv(cg,A);
2208: l=(long)((pariC2*(prec-2) + fabs(gtodouble(s)-1.)*log(rtodbl(cga)))
2209: / rtodbl(cgb)+1);
2210: v = anell(e, min(l,TEMPMAX));
2211: if (!flun) { s2=gsubsg(2,s); ns=gpui(cg,gsubgs(gmul2n(s,1),2),prec); }
2212: z=gzero;
2213: if (typ(s)==t_INT)
2214: {
2215: if (signe(s)<=0) { avma=av; return gzero; }
2216: gs=mpfactr(itos(s)-1,prec);
2217: }
2218: else gs=ggamma(s,prec);
2219: av1=avma; lim=stack_lim(av1,1);
2220: for (n=1; n<=l; n++)
2221: {
2222: p1=gdiv(incgam4(s,gmulsg(n,cga),gs,prec),gpui(stoi(n),s,prec));
2223: p2=flun? p1: gdiv(gmul(ns,incgam(s2,gmulsg(n,cgb),prec)),
2224: gpui(stoi(n),s2,prec));
2225: if (eps<0) p2=gneg_i(p2);
2226: z = gadd(z,gmul(gadd(p1,p2),(n<=TEMPMAX)? (GEN)v[n]: akell(e,stoi(n))));
2227: if (low_stack(lim, stack_lim(av1,1)))
2228: {
2229: if(DEBUGMEM>1) err(warnmem,"lseriesell");
2230: tetpil=avma; z=gerepile(av1,tetpil,gcopy(z));
2231: }
2232: }
2233: tetpil=avma; return gerepile(av,tetpil,gdiv(z,gs));
2234: }
2235:
2236: /********************************************************************/
2237: /** **/
2238: /** Tate's algorithm e (cf Anvers IV) **/
2239: /** Kodaira types, global minimal model **/
2240: /** **/
2241: /********************************************************************/
2242:
2243: /* Given an integral elliptic curve in ellinit form, and a prime p, returns the
2244: type of the fiber at p of the Neron model, as well as the change of variables
2245: in the form [f, kod, v, c].
2246:
2247: * The integer f is the conductor's exponent.
2248:
2249: * The integer kod is the Kodaira type using the following notation:
2250: II , III , IV --> 2, 3, 4
2251: I0 --> 1
2252: Inu --> 4+nu for nu > 0
2253: A '*' negates the code (e.g I* --> -2)
2254:
2255: * v is a quadruple [u, r, s, t] yielding a minimal model
2256:
2257: * c is the Tamagawa number.
2258:
2259: Uses Tate's algorithm (Anvers IV). Given the remarks at the bottom of
2260: page 46, the "long" algorithm is used for p < 4 only. */
2261: static void cumule(GEN *vtotal, GEN *e, GEN u, GEN r, GEN s, GEN t);
2262: static void cumule1(GEN *vtotal, GEN *e, GEN v2);
2263:
2264: static GEN
2265: localreduction_result(long av, long f, long kod, long c, GEN v)
2266: {
2267: long tetpil = avma;
2268: GEN result = cgetg(5, t_VEC);
2269: result[1] = lstoi(f); result[2] = lstoi(kod);
2270: result[3] = lcopy(v); result[4] = lstoi(c);
2271: return gerepile(av,tetpil, result);
2272: }
2273:
2274: /* ici, p1 != 2 et p1 != 3 */
2275: static GEN
2276: localreduction_carac_not23(GEN e, GEN p1)
2277: {
2278: long av = avma, k, f, kod, c, nuj, nudelta;
2279: GEN pk, p2k, a2prime, a3prime;
2280: GEN p2, r = gzero, s = gzero, t = gzero, v;
2281: GEN c4, c6, delta, unmodp, xun, tri, var, p4k, p6k;
2282:
2283: nudelta = ggval((GEN)e[12], p1);
2284: v = cgetg(5,t_VEC); v[1] = un; v[2] = v[3] = v[4] = zero;
2285: nuj = gcmp0((GEN)e[13]) ? 0 : - ggval((GEN)e[13], p1);
2286: k = (nuj > 0 ? nudelta - nuj : nudelta) / 12;
2287: c4 = (GEN)e[10]; c6 = (GEN)e[11]; delta = (GEN)e[12];
2288: if (k > 0) /* modele non minimal */
2289: {
2290: pk = gpuigs(p1, k);
2291: if (mpodd((GEN)e[1]))
2292: s = shifti(subii(pk, (GEN)e[1]), -1);
2293: else
2294: s = negi(shifti((GEN)e[1], -1));
2295: p2k = sqri(pk);
2296: p4k = sqri(p2k);
2297: p6k = mulii(p4k, p2k);
2298:
2299: a2prime = subii((GEN)e[2], mulii(s, addii((GEN)e[1], s)));
2300: switch(smodis(a2prime, 3))
2301: {
2302: case 0: r = negi(divis(a2prime, 3)); break;
2303: case 1: r = divis(subii(p2k, a2prime), 3); break;
2304: case 2: r = negi(divis(addii(a2prime, p2k), 3)); break;
2305: }
2306: a3prime = ellLHS0_i(e,r);
2307: if (mpodd(a3prime))
2308: t = shifti(subii(mulii(pk, p2k), a3prime), -1);
2309: else
2310: t = negi(shifti(a3prime, -1));
2311: v[1] = (long)pk; v[2] = (long)r; v[3] = (long)s; v[4] = (long)t;
2312: nudelta -= 12 * k;
2313: c4 = divii(c4, p4k); c6 = divii(c6, p6k);
2314: delta = divii(delta, sqri(p6k));
2315: }
2316: if (nuj > 0) switch(nudelta - nuj)
2317: {
2318: case 0: f = 1; kod = 4+nuj; /* Inu */
2319: switch(kronecker(negi(c6),p1))
2320: {
2321: case 1: c = nudelta; break;
2322: case -1: c = 2 - (nudelta % 2); break;
2323: default: err(tater1);
2324: }
2325: break;
2326: case 6: f = 2; kod = -4-nuj; /* Inu* */
2327: if (nuj & 1)
2328: c = 3 + kronecker(divii(mulii(c6, delta),gpuigs(p1, 9+nuj)), p1);
2329: else
2330: c = 3 + kronecker(divii(delta, gpuigs(p1, 6+nuj)), p1);
2331: break;
2332: default: err(tater1);
2333: }
2334: else switch(nudelta)
2335: {
2336: case 0: f = 0; kod = 1; c = 1; break; /* I0, regulier */
2337: case 2: f = 2; kod = 2; c = 1; break; /* II */
2338: case 3: f = 2; kod = 3; c = 2; break; /* III */
2339: case 4: f = 2; kod = 4; /* IV */
2340: c = 2 + kronecker(gdiv(mulis(c6, -6), sqri(p1)), p1);
2341: break;
2342: case 6: f = 2; kod = -1; /* I0* */
2343: p2 = sqri(p1);
2344: unmodp = gmodulsg(1,p1);
2345: var = gmul(unmodp,polx[0]);
2346: tri = gsub(gsqr(var),gmul(divii(gmulsg(3, c4), p2),unmodp));
2347: tri = gsub(gmul(tri, var),
2348: gmul(divii(gmul2n(c6,1), mulii(p2,p1)),unmodp));
2349: xun = gmodulcp(var,tri);
2350: c = lgef(ggcd((GEN)(gsub(gpui(xun,p1,0),xun))[2], tri)) - 2;
2351: break;
2352: case 8: f = 2; kod = -4; /* IV* */
2353: c = 2 + kronecker(gdiv(mulis(c6,-6), gpuigs(p1,4)), p1);
2354: break;
2355: case 9: f = 2; kod = -3; c = 2; break; /* III* */
2356: case 10: f = 2; kod = -2; c = 1; break; /* II* */
2357: default: err(tater1);
2358: }
2359: return localreduction_result(av,f,kod,c,v);
2360: }
2361:
2362: /* renvoie a_{ k,l } avec les notations de Tate */
2363: static int
2364: aux(GEN ak, int p, int l)
2365: {
2366: long av = avma, pl = p, res;
2367: while (--l) pl *= p;
2368: res = smodis(divis(ak, pl), p);
2369: avma = av; return res;
2370: }
2371:
2372: static int
2373: aux2(GEN ak, int p, GEN pl)
2374: {
2375: long av = avma, res;
2376: res = smodis(divii(ak, pl), p);
2377: avma = av;
2378: return res;
2379: }
2380:
2381: /* renvoie le nombre de racines distinctes du polynome XXX + aXX + bX + c
2382: * modulo p s'il y a une racine multiple, elle est renvoyee dans *mult
2383: */
2384: static int
2385: numroots3(int a, int b, int c, int p, int *mult)
2386: {
2387: if (p == 2)
2388: {
2389: if ((c + a * b) & 1) return 3;
2390: else { *mult = b; return (a + b) & 1 ? 2 : 1; }
2391: }
2392: else
2393: {
2394: if (a % 3) { *mult = a * b; return (a * b * (1 - b) + c) % 3 ? 3 : 2; }
2395: else { *mult = -c; return b % 3 ? 3 : 1; }
2396: }
2397: }
2398:
2399: /* idem pour aXX +bX + c */
2400: static int
2401: numroots2(int a, int b, int c, int p, int *mult)
2402: {
2403: if (p == 2) { *mult = c; return b & 1 ? 2 : 1; }
2404: else { *mult = a * b; return (b * b - a * c) % 3 ? 2 : 1; }
2405: }
2406:
2407: /* ici, p1 = 2 ou p1 = 3 */
2408: static GEN
2409: localreduction_carac_23(GEN e, GEN p1)
2410: {
2411: long av = avma, p, c, nu, nudelta;
2412: int a21, a42, a63, a32, a64, theroot, al, be, ga;
2413: GEN pk, p2k, pk1, p4, p6;
2414: GEN p2, p3, r = gzero, s = gzero, t = gzero, v;
2415:
2416: nudelta = ggval((GEN)e[12], p1);
2417: v = cgetg(5,t_VEC); v[1] = un; v[2] = v[3] = v[4] = zero;
2418:
2419: for(;;)
2420: {
2421: if (!nudelta)
2422: return localreduction_result(av, 0, 1, 1, v);
2423: /* I0 */
2424: p = itos(p1);
2425: if (!divise((GEN)e[6], p1))
2426: {
2427: if (smodis(negi((GEN)e[11]), p == 2 ? 8 : 3) == 1)
2428: c = nudelta;
2429: else
2430: c = 2 - (nudelta & 1);
2431: return localreduction_result(av, 1, 4 + nudelta, c, v);
2432: }
2433: /* Inu */
2434: if (p == 2)
2435: {
2436: r = modis((GEN)e[4], 2);
2437: s = modis(addii(r, (GEN)e[2]), 2);
2438: if (signe(r)) t = modis(addii(addii((GEN)e[4], (GEN)e[5]), s), 2);
2439: else t = modis((GEN)e[5], 2);
2440: }
2441: else /* p == 3 */
2442: {
2443: r = negi(modis((GEN)e[8], 3));
2444: s = modis((GEN)e[1], 3);
2445: t = modis(ellLHS0_i(e,r), 3);
2446: }
2447: cumule(&v, &e, gun, r, s, t); /* p | a1, a2, a3, a4 et a6 */
2448: p2 = stoi(p*p);
2449: if (!divise((GEN)e[5], p2))
2450: return localreduction_result(av, nudelta, 2, 1, v);
2451: /* II */
2452: p3 = stoi(p*p*p);
2453: if (!divise((GEN)e[9], p3))
2454: return localreduction_result(av, nudelta - 1, 3, 2, v);
2455: /* III */
2456: if (!divise((GEN)e[8], p3))
2457: {
2458: if (smodis((GEN)e[8], (p==2)? 32: 27) == p*p)
2459: c = 3;
2460: else
2461: c = 1;
2462: return localreduction_result(av, nudelta - 2, 4, c, v);
2463: }
2464: /* IV */
2465:
2466: /* now for the last five cases... */
2467:
2468: if (!divise((GEN)e[5], p3))
2469: cumule(&v, &e, gun, gzero, gzero, p == 2? gdeux: modis((GEN)e[3], 9));
2470: /* p | a1, a2; p^2 | a3, a4; p^3 | a6 */
2471: a21 = aux((GEN)e[2], p, 1); a42 = aux((GEN)e[4], p, 2);
2472: a63 = aux((GEN)e[5], p, 3);
2473: switch (numroots3(a21, a42, a63, p, &theroot))
2474: {
2475: case 3:
2476: if (p == 2)
2477: c = 1 + (a63 == 0) + ((a21 + a42 + a63) & 1);
2478: else
2479: c = 1 + (a63 == 0) + (((1 + a21 + a42 + a63) % 3) == 0)
2480: + (((1 - a21 + a42 - a63) % 3) == 0);
2481: return localreduction_result(av, nudelta - 4, -1, c, v);
2482: /* I0* */
2483: case 2: /* calcul de nu */
2484: if (theroot) cumule(&v, &e, gun, stoi(theroot * p), gzero, gzero);
2485: /* p | a1; p^2 | a2, a3; p^3 | a4; p^4 | a6 */
2486: nu = 1;
2487: pk = p2;
2488: p2k = stoi(p * p * p * p);
2489: for(;;)
2490: {
2491: if (numroots2(al = 1, be = aux2((GEN)e[3], p, pk),
2492: ga = -aux2((GEN)e[5], p, p2k), p, &theroot) == 2)
2493: break;
2494: if (theroot) cumule(&v, &e, gun, gzero, gzero, mulsi(theroot,pk));
2495: pk1 = pk; pk = mulsi(p, pk); p2k = mulsi(p, p2k);
2496: nu++;
2497: if (numroots2(al = a21, be = aux2((GEN)e[4], p, pk),
2498: ga = aux2((GEN)e[5], p, p2k), p, &theroot) == 2)
2499: break;
2500: if (theroot) cumule(&v, &e, gun, mulsi(theroot, pk1), gzero, gzero);
2501: p2k = mulsi(p, p2k);
2502: nu++;
2503: }
2504: if (p == 2)
2505: c = 4 - 2 * (ga & 1);
2506: else
2507: c = 3 + kross(be * be - al * ga, 3);
2508: return localreduction_result(av, nudelta - 4 - nu, -4 - nu, c, v);
2509: /* Inu* */
2510: case 1:
2511: if (theroot) cumule(&v, &e, gun, stoi(theroot * p), gzero, gzero);
2512: /* p | a1; p^2 | a2, a3; p^3 | a4; p^4 | a6 */
2513: a32 = aux((GEN)e[3], p, 2); a64 = aux((GEN)e[5], p, 4);
2514: if (numroots2(1, a32, -a64, p, &theroot) == 2)
2515: {
2516: if (p == 2)
2517: c = 3 - 2 * a64;
2518: else
2519: c = 2 + kross(a32 * a32 + a64, 3);
2520: return localreduction_result(av, nudelta - 6, -4, c, v);
2521: }
2522: /* IV* */
2523: if (theroot) cumule(&v, &e, gun, gzero, gzero, stoi(theroot*p*p));
2524: /* p | a1; p^2 | a2; p^3 | a3, a4; p^5 | a6 */
2525: p4 = sqri(p2);
2526: if (!divise((GEN)e[4], p4))
2527: return localreduction_result(av, nudelta - 7, -3, 2, v);
2528: /* III* */
2529: p6 = mulii(p4, p2);
2530: if (!divise((GEN)e[5], p6))
2531: return localreduction_result(av, nudelta - 8, -2, 1, v);
2532: /* II* */
2533: cumule(&v, &e, p1, gzero, gzero, gzero); /* non minimal, on repart
2534: pour un tour */
2535: nudelta -= 12;
2536: }
2537: }
2538: /* Not reached */
2539: }
2540:
2541: GEN
2542: localreduction(GEN e, GEN p1)
2543: {
2544: checkell(e);
2545: if (typ(e[12]) != t_INT)
2546: err(talker,"not an integral curve in localreduction");
2547: if (gcmpgs(p1, 3) > 0) /* p different de 2 ou 3 */
2548: return localreduction_carac_not23(e,p1);
2549: else
2550: return localreduction_carac_23(e,p1);
2551: }
2552:
2553: #if 0
2554: /* Calcul de toutes les fibres non elliptiques d'une courbe sur Z.
2555: * Etant donne une courbe elliptique sous forme longue e, dont les coefficients
2556: * sont entiers, renvoie une matrice dont les lignes sont de la forme
2557: * [p, fp, kodp, cp]. Il y a une ligne par diviseur premier du discriminant.
2558: */
2559: GEN
2560: globaltatealgo(GEN e)
2561: {
2562: long k, l,av;
2563: GEN p1, p2, p3, p4, prims, result;
2564:
2565: checkell(e);
2566: prims = decomp((GEN)e[12]);
2567: l = lg(p1 = (GEN)prims[1]);
2568: p2 = (GEN)prims[2];
2569: if ((long)prims == avma) cgiv(prims);
2570: result = cgetg(5, t_MAT);
2571: result[1] = (long)p1;
2572: result[2] = (long)p2;
2573: result[3] = (long)(p3 = cgetg(l, t_COL));
2574: for (k = 1; k < l; k++) p3[k] = lgeti(3);
2575: result[4] = (long)(p4 = cgetg(l, t_COL));
2576: for (k = 1; k < l; k++) p4[k] = lgeti(3);
2577: av = avma;
2578: for (k = 1; k < l; k++)
2579: {
2580: GEN q = localreduction(e, (GEN)p1[k]);
2581: affii((GEN)q[1],(GEN)p2[k]);
2582: affii((GEN)q[2],(GEN)p3[k]);
2583: affii((GEN)q[4],(GEN)p4[k]);
2584: avma = av;
2585: }
2586: return result;
2587: }
2588: #endif
2589:
2590: /* Algorithme de reduction d'une courbe sur Q a sa forme standard. Etant
2591: * donne une courbe elliptique sous forme longue e, dont les coefficients
2592: * sont rationnels, renvoie son [N, [u, r, s, t], c], ou N est le conducteur
2593: * arithmetique de e, [u, r, s, t] est le changement de variables qui reduit
2594: * e a sa forme minimale globale dans laquelle a1 et a3 valent 0 ou 1, et a2
2595: * vaut -1, 0 ou 1 et tel que u est un rationnel positif. Enfin c est le
2596: * produit des nombres de Tamagawa locaux cp.
2597: */
2598: GEN
2599: globalreduction(GEN e1)
2600: {
2601: long i, k, l, m, tetpil, av = avma;
2602: GEN p1, c = gun, prims, result, N = gun, u = gun, r, s, t;
2603: GEN v = cgetg(5, t_VEC), a = cgetg(7, t_VEC), e = cgetg(20, t_VEC);
2604:
2605: checkell(e1);
2606: for (i = 1; i < 5; i++) a[i] = e1[i]; a[5] = zero; a[6] = e1[5];
2607: prims = decomp(denom(a));
2608: p1 = (GEN)prims[1]; l = lg(p1);
2609: for (k = 1; k < l; k++)
2610: {
2611: int n = 0;
2612: for (i = 1; i < 7; i++)
2613: if (!gcmp0((GEN)a[i]))
2614: {
2615: m = i * n + ggval((GEN)a[i], (GEN)p1[k]);
2616: while (m < 0) { n++; m += i; }
2617: }
2618: u = gmul(u, gpuigs((GEN)p1[k], n));
2619: }
2620: v[1] = linv(u); v[2] = v[3] = v[4] = zero;
2621: for (i = 1; i < 14; i++) e[i] = e1[i];
2622: for (; i < 20; i++) e[i] = zero;
2623: if (!gcmp1(u)) e = coordch(e, v);
2624: prims = decomp((GEN)e[12]);
2625: l = lg(p1 = (GEN)prims[1]);
2626: for (k = (signe(e[12]) < 0) + 1; k < l; k++)
2627: {
2628: GEN q = localreduction(e, (GEN)p1[k]);
2629: GEN v1 = (GEN)q[3];
2630: N = mulii(N, gpui((GEN)p1[k],(GEN)q[1],0));
2631: c = mulii(c, (GEN)q[4]);
2632: if (!gcmp1((GEN)v1[1])) cumule1(&v, &e, v1);
2633: }
2634: s = gdiventgs((GEN)e[1], -2);
2635: r = gdiventgs(gaddgs(gsub(gsub((GEN)e[2], gmul(s,(GEN)e[1])), gsqr(s)), 1), -3);
2636: t = gdiventgs(ellLHS0(e,r), -2);
2637: cumule(&v, &e, gun, r, s, t);
2638: tetpil = avma;
2639: result = cgetg(4, t_VEC); result[1] = lcopy(N); result[2] = lcopy(v);
2640: result[3] = lcopy(c);
2641: return gerepile(av, tetpil, result);
2642: }
2643:
2644: /* cumule les effets de plusieurs chgts de variable. On traite a part les cas
2645: * particuliers frequents, tels que soit u = 1, soit r' = s' = t' = 0
2646: */
2647: static void
2648: cumulev(GEN *vtotal, GEN u, GEN r, GEN s, GEN t)
2649: {
2650: long av = avma, tetpil;
2651: GEN temp, v = *vtotal, v3 = cgetg(5, t_VEC);
2652: if (gcmp1((GEN)v[1]))
2653: {
2654: v3[1] = lcopy(u);
2655: v3[2] = ladd((GEN)v[2], r);
2656: v3[3] = ladd((GEN)v[3], s);
2657: av = avma;
2658: temp = gadd((GEN)v[4], gmul((GEN)v[3], r));
2659: tetpil = avma;
2660: v3[4] = lpile(av, tetpil, gadd(temp, t));
2661: }
2662: else if (gcmp0(r) && gcmp0(s) && gcmp0(t))
2663: {
2664: v3[1] = lmul((GEN)v[1], u);
2665: v3[2] = lcopy((GEN)v[2]);
2666: v3[3] = lcopy((GEN)v[3]);
2667: v3[4] = lcopy((GEN)v[4]);
2668: }
2669: else /* cas general */
2670: {
2671: v3[1] = lmul((GEN)v[1], u);
2672: temp = gsqr((GEN)v[1]);
2673: v3[2] = ladd(gmul(temp, r), (GEN)v[2]);
2674: v3[3] = ladd(gmul((GEN)v[1], s), (GEN)v[3]);
2675: v3[4] = ladd((GEN)v[4], gmul(temp, gadd(gmul((GEN)v[1], t), gmul((GEN)v[3], r))));
2676:
2677: tetpil = avma;
2678: v3 = gerepile(av, tetpil, gcopy(v3));
2679: }
2680: *vtotal = v3;
2681: }
2682:
2683: static void
2684: cumule(GEN *vtotal, GEN *e, GEN u, GEN r, GEN s, GEN t)
2685: {
2686: long av = avma, tetpil;
2687: GEN v2 = cgetg(5, t_VEC);
2688: v2[1] = (long)u; v2[2] = (long)r; v2[3] = (long)s; v2[4] = (long)t;
2689: tetpil = avma;
2690: *e = gerepile(av, tetpil, coordch(*e, v2));
2691: cumulev(vtotal, u, r, s, t);
2692: }
2693:
2694: static void
2695: cumule1(GEN *vtotal, GEN *e, GEN v2)
2696: {
2697: *e = coordch(*e, v2);
2698: cumulev(vtotal, (GEN)v2[1], (GEN)v2[2], (GEN)v2[3], (GEN)v2[4]);
2699: }
2700:
2701: /********************************************************************/
2702: /** **/
2703: /** Parametrisation modulaire **/
2704: /** **/
2705: /********************************************************************/
2706:
2707: GEN
2708: taniyama(GEN e)
2709: {
2710: GEN v,w,c,d,s1,s2,s3;
2711: long n,m,av=avma,tetpil;
2712:
2713: checkell(e); v = cgetg(precdl+3,t_SER);
2714: v[1] = evalsigne(1) | evalvalp(-2) | evalvarn(0);
2715: v[2] = un;
2716: c=gtoser(anell(e,precdl+1),0); setvalp(c,1);
2717: d=ginv(c); c=gsqr(d);
2718: for (n=-3; n<=precdl-4; n++)
2719: {
2720: if (n!=2)
2721: {
2722: s3=n?gzero:(GEN)e[7];
2723: if (n>-3) s3=gadd(s3,gmul((GEN)e[6],(GEN)v[n+4]));
2724: s2=gzero;
2725: for (m=-2; m<=n+1; m++)
2726: s2 = gadd(s2,gmulsg(m*(n+m),gmul((GEN)v[m+4],(GEN)c[n-m+4])));
2727: s2=gmul2n(s2,-1);
2728: s1=gzero;
2729: for (m=-1; m+m<=n; m++)
2730: {
2731: if (m+m==n)
2732: s1=gadd(s1,gsqr((GEN)v[m+4]));
2733: else
2734: s1=gadd(s1,gmul2n(gmul((GEN)v[m+4],(GEN)v[n-m+4]),1));
2735: }
2736: v[n+6]=ldivgs(gsub(gadd(gmulsg(6,s1),s3),s2),(n+2)*(n+1)-12);
2737: }
2738: else
2739: {
2740: setlg(v,9); v[8]=(long)polx[MAXVARN];
2741: w=deriv(v,0); setvalp(w,-2);
2742: s1=gadd((GEN)e[8],gmul(v,gadd(gmul2n((GEN)e[7],1),gmul(v,gadd((GEN)e[6],gmul2n(v,2))))));
2743: setlg(v,precdl+3);
2744: s2=gsub(s1,gmul(c,gsqr(w)));
2745: s2=gsubst((GEN)s2[2],MAXVARN,polx[0]);
2746: v[n+6]=lneg(gdiv((GEN)s2[2],(GEN)s2[3]));
2747: }
2748: }
2749: w=gsub(gmul(polx[0],gmul(d,deriv(v,0))), ellLHS0(e,v));
2750: tetpil=avma; s1=cgetg(3,t_VEC); s1[1]=lcopy(v); s1[2]=lmul2n(w,-1);
2751: return gerepile(av,tetpil,s1);
2752: }
2753:
2754: /********************************************************************/
2755: /** **/
2756: /** TORSION POINTS (over Q) **/
2757: /** **/
2758: /********************************************************************/
2759: /* assume e is defined over Q (use Mazur's theorem) */
2760: GEN
2761: orderell(GEN e, GEN p)
2762: {
2763: GEN p1;
2764: long av=avma,k;
2765:
2766: checkell(e); checkpt(p);
2767: k=typ(e[13]);
2768: if (k!=t_INT && !is_frac_t(k))
2769: err(impl,"orderell for nonrational elliptic curves");
2770: p1=p; k=1;
2771: for (k=1; k<16; k++)
2772: {
2773: if (lg(p1)<3) { avma=av; return stoi(k); }
2774: p1 = addell(e,p1,p);
2775: }
2776: avma=av; return gzero;
2777: }
2778:
2779: /* one can do much better by factoring denom(D) (see ellglobalred) */
2780: static GEN
2781: ellintegralmodel(GEN e)
2782: {
2783: GEN a = cgetg(6,t_VEC), v;
2784: long i;
2785:
2786: for (i=1; i<6; i++) a[i]=e[i];
2787: a = denom(a); if (gcmp1(a)) return NULL;
2788: v = cgetg(5,t_VEC);
2789: v[1]=linv(a); v[2]=v[3]=v[4]=zero; return v;
2790: }
2791:
2792: /* Using Lutz-Nagell */
2793:
2794: /* p is a polynomial of degree exactly 3 with integral coefficients
2795: * and leading term 4. Outputs the vector of rational roots of p
2796: */
2797: static GEN
2798: ratroot(GEN p)
2799: {
2800: GEN v,a,ld;
2801: long i,t;
2802:
2803: i=2; while (!signe(p[i])) i++;
2804: if (i==5)
2805: { v=cgetg(2,t_VEC); v[1]=zero; return v; }
2806: if (i==4)
2807: { v=cgetg(3,t_VEC); v[1]=zero; v[2]=ldivgs((GEN)p[4],-4); return v; }
2808:
2809: v=cgetg(4,t_VEC); t=1;
2810: if (i==3) v[t++]=zero;
2811: ld=divisors(gmul2n((GEN)p[i],2));
2812: for (i=1; i<lg(ld); i++)
2813: {
2814: a = gmul2n((GEN)ld[i],-2);
2815: if (!gsigne(poleval(p,a))) v[t++]=(long)a;
2816: a = gneg_i(a);
2817: if (!gsigne(poleval(p,a))) v[t++]=(long)a;
2818: }
2819: setlg(v,t); return v;
2820: }
2821:
2822: static int
2823: is_new_torsion(GEN e, GEN v, GEN p, long t2) {
2824: GEN pk = p, pkprec = NULL;
2825: long k,l;
2826:
2827: for (k=2; k<=6; k++)
2828: {
2829: pk=addell(e,pk,p);
2830: if (lg(pk)==2) return 1;
2831:
2832: for (l=2; l<=t2; l++)
2833: if (gegal((GEN)pk[1],gmael(v,l,1))) return 1;
2834:
2835: if (pkprec && k<=5)
2836: if (gegal((GEN)pk[1],(GEN)pkprec[1])) return 1;
2837: pkprec=pk;
2838: }
2839: return 0;
2840: }
2841:
2842: GEN
2843: torsellnagelllutz(GEN e)
2844: {
2845: GEN d,ld,pol,p1,lr,r,v,w,w2,w3;
2846: long i,j,nlr,t,t2,k,k2,av=avma;
2847:
2848: checkell(e);
2849: v = ellintegralmodel(e);
2850: if (v) e = coordch(e,v);
2851: pol = RHSpol(e);
2852: lr=ratroot(pol); nlr=lg(lr)-1;
2853: r=cgetg(17,t_VEC); p1=cgetg(2,t_VEC); p1[1]=zero; r[1]=(long)p1;
2854: for (t=1,i=1; i<=nlr; i++)
2855: {
2856: p1=cgetg(3,t_VEC);
2857: p1[1] = lr[i];
2858: p1[2] = lmul2n(gneg(ellLHS0(e,(GEN)lr[i])), -1);
2859: r[++t]=(long)p1;
2860: }
2861: ld = factor(gmul2n(absi((GEN)e[12]), 4));
2862: p1 = (GEN)ld[2]; k = lg(p1);
2863: for (i=1; i<k; i++) p1[i] = lshifti((GEN)p1[i], -1);
2864: ld = divisors(ld);
2865: for (t2=t,j=1; j<lg(ld); j++)
2866: {
2867: d=(GEN)ld[j]; lr=ratroot(gsub(pol,gsqr(d)));
2868: for (i=1; i<lg(lr); i++)
2869: {
2870: p1 = cgetg(3,t_VEC);
2871: p1[1] = lr[i];
2872: p1[2] = lmul2n(gsub(d,ellLHS0(e,(GEN)lr[i])), -1);
2873:
2874: if (is_new_torsion(e,r,p1,t2))
2875: {
2876: GEN p2 = cgetg(3,t_VEC);
2877: p2[1] = p1[1];
2878: p2[2] = lsub((GEN)p1[2],d);
2879: r[++t]=(long)p1;
2880: r[++t]=(long)p2;
2881: }
2882: }
2883: }
2884: if (t==1)
2885: {
2886: avma=av; w=cgetg(4,t_VEC);
2887: w[1] = un;
2888: w[2] = lgetg(1,t_VEC);
2889: w[3] = lgetg(1,t_VEC);
2890: return w;
2891: }
2892:
2893: if (nlr<3)
2894: {
2895: w2=cgetg(2,t_VEC); w2[1]=lstoi(t);
2896: for (k=2; k<=t; k++)
2897: if (itos(orderell(e,(GEN)r[k])) == t) break;
2898: if (k>t) err(bugparier,"torsell (bug1)");
2899:
2900: w3=cgetg(2,t_VEC); w3[1]=r[k];
2901: }
2902: else
2903: {
2904: if (t&3) err(bugparier,"torsell (bug2)");
2905: t2 = t>>1;
2906: w2=cgetg(3,t_VEC); w2[1]=lstoi(t2); w2[2]=(long)gdeux;
2907: for (k=2; k<=t; k++)
2908: if (itos(orderell(e,(GEN)r[k])) == t2) break;
2909: if (k>t) err(bugparier,"torsell (bug3)");
2910:
2911: p1 = powell(e,(GEN)r[k],stoi(t>>2));
2912: k2 = (lg(p1)==3 && gegal((GEN)r[2],p1))? 3: 2;
2913: w3=cgetg(3,t_VEC); w3[1]=r[k]; w3[2]=r[k2];
2914: }
2915: if (v)
2916: {
2917: v[1] = linv((GEN)v[1]);
2918: w3 = pointch(w3,v);
2919: }
2920: w=cgetg(4,t_VEC);
2921: w[1] = lstoi(t);
2922: w[2] = (long)w2;
2923: w[3] = (long)w3;
2924: return gerepileupto(av, gcopy(w));
2925: }
2926:
2927: /* Using Doud's algorithm */
2928:
2929: /* Input e and n, finds a bound for #Tor */
2930: static long
2931: torsbound(GEN e, long n)
2932: {
2933: long av = avma, m, b, c, d, prime = 2;
2934: byteptr p = diffptr;
2935: GEN D = (GEN)e[12];
2936:
2937: b = c = m = 0; p++;
2938: while (m<n)
2939: {
2940: d = *p++; if (!d) err(primer1);
2941: prime += d;
2942: if (ggval(D,stoi(prime)) == 0)
2943: {
2944: b = cgcd(b, prime+1 - itos(apell0(e,prime)));
2945: if (b==c) m++; else {c = b; m = 0;}
2946: avma = av;
2947: }
2948: }
2949: return b;
2950: }
2951:
2952: static GEN
2953: _round(GEN x, long *e)
2954: {
2955: GEN y = grndtoi(x,e);
2956: if (*e > -5 && bit_accuracy(gprecision(x)) < gexpo(y) - 10)
2957: err(talker, "ellinit data not accurate enough. Increase precision");
2958: return y;
2959: }
2960:
2961: /* Input the curve, a point, and an integer n, returns a point of order n
2962: on the curve, or NULL if q is not rational. */
2963: static GEN
2964: torspnt(GEN E, GEN q, long n)
2965: {
2966: GEN p = cgetg(3,t_VEC);
2967: long e;
2968: p[1] = lmul2n(_round(gmul2n((GEN)q[1],2), &e),-2);
2969: if (e > -5) return NULL;
2970: p[2] = lmul2n(_round(gmul2n((GEN)q[2],3), &e),-3);
2971: if (e > -5) return NULL;
2972: return (gcmp0(gimag(p)) && oncurve(E,p)
2973: && lg(powell(E,p,stoi(n))) == 2
2974: && itos(orderell(E,p)) == n)? greal(p): NULL;
2975: }
2976:
2977: static int
2978: smaller_x(GEN p, GEN q)
2979: {
2980: int s = absi_cmp(denom(p), denom(q));
2981: return (s<0 || (s==0 && absi_cmp(numer(p),numer(q)) < 0));
2982: }
2983:
2984: /* best generator in cycle of length k */
2985: static GEN
2986: best_in_cycle(GEN e, GEN p, long k)
2987: {
2988: GEN p0 = p,q = p;
2989: long i;
2990:
2991: for (i=2; i+i<k; i++)
2992: {
2993: q = addell(e,q,p0);
2994: if (cgcd(i,k)==1 && smaller_x((GEN)q[1], (GEN)p[1])) p = q;
2995: }
2996: return (gsigne(d_ellLHS(e,p)) < 0)? invell(e,p): p;
2997: }
2998:
2999: static GEN
3000: tors(GEN e, long k, GEN p, GEN q, GEN v)
3001: {
3002: GEN p1,r;
3003: if (q)
3004: {
3005: long n = k>>1;
3006: GEN p1, best = q, np = powell(e,p,stoi(n));
3007: if (n % 2 && smaller_x((GEN)np[1], (GEN)best[1])) best = np;
3008: p1 = addell(e,q,np);
3009: if (smaller_x((GEN)p1[1], (GEN)best[1])) q = p1;
3010: else if (best == np) { p = addell(e,p,q); q = np; }
3011: p = best_in_cycle(e,p,k);
3012: if (v)
3013: {
3014: v[1] = linv((GEN)v[1]);
3015: p = pointch(p,v);
3016: q = pointch(q,v);
3017: }
3018: r = cgetg(4,t_VEC);
3019: r[1] = lstoi(2*k); p1 = cgetg(3,t_VEC); p1[1] = lstoi(k); p1[2] = deux;
3020: r[2] = (long)p1; p1 = cgetg(3,t_VEC); p1[1] = lcopy(p); p1[2] = lcopy(q);
3021: r[3] = (long)p1;
3022: }
3023: else
3024: {
3025: if (p)
3026: {
3027: p = best_in_cycle(e,p,k);
3028: if (v)
3029: {
3030: v[1] = linv((GEN)v[1]);
3031: p = pointch(p,v);
3032: }
3033: r = cgetg(4,t_VEC);
3034: r[1] = lstoi(k); p1 = cgetg(2,t_VEC); p1[1] = r[1];
3035: r[2] = (long)p1; p1 = cgetg(2,t_VEC); p1[1] = lcopy(p);
3036: r[3] = (long)p1;
3037: }
3038: else
3039: {
3040: r = cgetg(4,t_VEC);
3041: r[1] = un;
3042: r[2] = lgetg(1,t_VEC);
3043: r[3] = lgetg(1,t_VEC);
3044: }
3045: }
3046: return r;
3047: }
3048:
3049: GEN
3050: torselldoud(GEN e)
3051: {
3052: long b,i,ord,av=avma,prec, k = 1;
3053: GEN v,w1,w22,w1j,w12,p,tor1,tor2;
3054:
3055: checkbell(e);
3056: v = ellintegralmodel(e);
3057: if (v) e = coordch(e,v);
3058:
3059: prec=precision((GEN)e[15]);
3060: prec=max(prec,MEDDEFAULTPREC);
3061: b = torsbound(e,3);
3062: if (b==1) { avma=av; return tors(e,1,NULL,NULL, v); }
3063: w22 = gmul2n((GEN)e[16],-1);
3064: w1 = (GEN)e[15];
3065: if (b % 4)
3066: {
3067: for (i=10; i>1; i--)
3068: {
3069: if (b%i==0)
3070: {
3071: w1j = gdivgs(w1,i);
3072: p = torspnt(e,pointell(e,w1j,prec),i);
3073: if (i%2==0 && p==NULL)
3074: {
3075: p = torspnt(e,pointell(e,gadd(w22,w1j),prec),i);
3076: if (p==NULL)
3077: p = torspnt(e,pointell(e,gadd(w22,gmul2n(w1j,1)),prec),i);
3078: }
3079: }
3080: else p = NULL;
3081: if (p) {k = i; break; }
3082: }
3083: return gerepileupto(av, tors(e,k,p,NULL, v));
3084: }
3085:
3086: ord = 0; tor2 = NULL;
3087: w12 = gmul2n((GEN)e[15],-1);
3088: if ((p = torspnt(e,pointell(e,w12,prec),2)))
3089: {
3090: tor1 = p; ord++;
3091: }
3092: if ((p = torspnt(e,pointell(e,w22,prec),2))
3093: || (!ord && (p = torspnt(e,pointell(e,gadd(w12,w22),prec),2))))
3094: {
3095: tor2 = p; ord += 2;
3096: }
3097:
3098: switch(ord)
3099: {
3100: case 0:
3101: for (i=9; i>1; i-=2)
3102: {
3103: w1j=gdivgs((GEN)e[15],i);
3104: p = (b%i==0)? torspnt(e,pointell(e,w1j,prec),i): NULL;
3105: if (p) { k = i; break; }
3106: }
3107: break;
3108:
3109: case 1:
3110: for (i=12; i>2; i-=2)
3111: {
3112: w1j=gdivgs((GEN)e[15],i);
3113: if (b%i==0)
3114: {
3115: p = torspnt(e,pointell(e,w1j,prec),i);
3116: if (i%4==0 && p==NULL)
3117: p = torspnt(e,pointell(e,gadd(w22,w1j),prec),i);
3118: }
3119: else p = NULL;
3120: if (p) { k = i; break; }
3121: }
3122: if (!p) { p = tor1; k = 2; }
3123: break;
3124:
3125: case 2:
3126: for (i=5; i>1; i-=2)
3127: {
3128: w1j = gdivgs((GEN)e[15],i);
3129: p = (b%i==0)? torspnt(e,pointell(e,gadd(w22,w1j),prec),i+i): NULL;
3130: if (p) { k = i+i; break; }
3131: }
3132: if (!p) { p = tor2; k = 2; }
3133: tor2 = NULL; break;
3134:
3135: case 3:
3136: for (i=8; i>2; i-=2)
3137: {
3138: w1j=gdivgs((GEN)e[15],i);
3139: p = (b%(2*i)==0)? torspnt(e,pointell(e,w1j,prec),i): NULL;
3140: if (p) { k = i; break; }
3141: }
3142: if (!p) { p = tor1; k = 2; }
3143: break;
3144: }
3145: return gerepileupto(av, tors(e,k,p,tor2, v));
3146: }
3147:
3148: GEN
3149: elltors0(GEN e, long flag)
3150: {
3151: switch(flag)
3152: {
3153: case 0: return torselldoud(e);
3154: case 1: return torsellnagelllutz(e);
3155: default: err(flagerr,"torsell");
3156: }
3157: return NULL; /* not reached */
3158: }
3159:
3160: /* par compatibilite */
3161: GEN torsell(GEN e) {return torselldoud(e);}
3162:
3163: /* LOCAL ROOT NUMBERS, D'APRES HALBERSTADT halberst@math.jussieu.fr */
3164:
3165: /* ici p=2 ou 3 */
3166: static long
3167: neron(GEN e, GEN p, long* ptkod)
3168: {
3169: long av=avma,kod,v4,v6,vd;
3170: GEN c4, c6, d, nv;
3171:
3172: nv=localreduction(e,p);
3173: kod=itos((GEN)nv[2]); *ptkod=kod;
3174: c4=(GEN)e[10]; c6=(GEN)e[11]; d=(GEN)e[12];
3175: v4=gcmp0(c4) ? 12 : ggval(c4,p);
3176: v6=gcmp0(c6) ? 12 : ggval(c6,p);
3177: vd=ggval(d,p);
3178: avma=av;
3179: switch(itos(p))
3180: {
3181: case 3:
3182: if (labs(kod)>4) return 1;
3183: else
3184: {
3185: switch(kod)
3186: {
3187: case -1: case 1: return v4&1 ? 2 : 1;
3188: case -3: case 3: return (2*v6>vd+3) ? 2 : 1;
3189: case -4: case 2:
3190: switch (vd%6)
3191: {
3192: case 4: return 3;
3193: case 5: return 4;
3194: default: return v6%3==1 ? 2 : 1;
3195: }
3196: default: /* kod = -2 et 4 */
3197: switch (vd%6)
3198: {
3199: case 0: return 2;
3200: case 1: return 3;
3201: default: return 1;
3202: }
3203: }
3204: }
3205: case 2:
3206: if (kod>4) return 1;
3207: else
3208: {
3209: switch(kod)
3210: {
3211: case 1: return (v6>0) ? 2 : 1;
3212: case 2:
3213: if (vd==4) return 1;
3214: else
3215: {
3216: if (vd==7) return 3;
3217: else return v4==4 ? 2 : 4;
3218: }
3219: case 3:
3220: switch(vd)
3221: {
3222: case 6: return 3;
3223: case 8: return 4;
3224: case 9: return 5;
3225: default: return v4==5 ? 2 : 1;
3226: }
3227: case 4: return v4>4 ? 2 : 1;
3228: case -1:
3229: switch(vd)
3230: {
3231: case 9: return 2;
3232: case 10: return 4;
3233: default: return v4>4 ? 3 : 1;
3234: }
3235: case -2:
3236: switch(vd)
3237: {
3238: case 12: return 2;
3239: case 14: return 3;
3240: default: return 1;
3241: }
3242: case -3:
3243: switch(vd)
3244: {
3245: case 12: return 2;
3246: case 14: return 3;
3247: case 15: return 4;
3248: default: return 1;
3249: }
3250: case -4: return v6==7 ? 2 : 1;
3251: case -5: return (v6==7 || v4==6) ? 2 : 1;
3252: case -6:
3253: switch(vd)
3254: {
3255: case 12: return 2;
3256: case 13: return 3;
3257: default: return v4==6 ? 2 : 1;
3258: }
3259: case -7: return (vd==12 || v4==6) ? 2 : 1;
3260: default: return v4==6 ? 2 : 1;
3261: }
3262: }
3263: default: return 0; /* should not occur */
3264: }
3265: }
3266:
3267: static long
3268: ellrootno_2(GEN e)
3269: {
3270: long n2,kod,u,v,x1,y1,d1,av=avma,v4,v6,w2;
3271: GEN p=gdeux,c4,c6,tmp,p6;
3272:
3273: n2=neron(e,p,&kod); c4=(GEN)e[10]; c6=(GEN)e[11]; p6=stoi(64);
3274: if (gcmp0(c4)) {v4=12; u=0;}
3275: else {v4=pvaluation(c4,p,&tmp); u=itos(modii(tmp,p6));}
3276: if (gcmp0(c6)) {v6=12; v=0;}
3277: else {v6=pvaluation(c6,p,&tmp); v=itos(modii(tmp,p6));}
3278: (void)pvaluation((GEN)e[12],p,&tmp); d1=itos(modii(tmp,p6));
3279: avma=av; x1=u+v+v;
3280: if (kod>=5)
3281: {w2=mpodd(addii((GEN)e[2],(GEN)e[3])) ? 1 : -1; avma=av; return w2;}
3282: if (kod<-9) return (n2==2) ? -kross(-1,v) : -1;
3283: switch(kod)
3284: {
3285: case 1: return 1;
3286: case 2:
3287: switch(n2)
3288: {
3289: case 1:
3290: switch(v4)
3291: {
3292: case 4: return kross(-1,u);
3293: case 5: return 1;
3294: default: return -1;
3295: }
3296: case 2: return (v6==7) ? 1 : -1;
3297: case 3: return (v%8==5 || (u*v)%8==5) ? 1 : -1;
3298: case 4: if (v4>5) return kross(-1,v);
3299: return (v4==5) ? -kross(-1,u) : -1;
3300: }
3301: case 3:
3302: switch(n2)
3303: {
3304: case 1: return -kross(2,u*v);
3305: case 2: return -kross(2,v);
3306: case 3: y1=itos(modis(gsubsg(u,gmul2n(c6,-5)),16)); avma=av;
3307: return (y1==7 || y1==11) ? 1 : -1;
3308: case 4: return (v%8==3 || (2*u+v)%8==7) ? 1 : -1;
3309: case 5: return v6==8 ? kross(2,x1) : kross(-2,u);
3310: }
3311: case -1:
3312: switch(n2)
3313: {
3314: case 1: return -kross(2,x1);
3315: case 2: return (v%8==7) || (x1%32==11) ? 1 : -1;
3316: case 3: return v4==6 ? 1 : -1;
3317: case 4: if (v4>6) return kross(-1,v);
3318: return v4==6 ? -kross(-1,u*v) : -1;
3319: }
3320: case -2: return n2==1 ? kross(-2,v) : kross(-1,v);
3321: case -3:
3322: switch(n2)
3323: {
3324: case 1: y1=(u-2*v)%64; if (y1<0) y1+=64;
3325: return (y1==3) || (y1==19) ? 1 : -1;
3326: case 2: return kross(2*kross(-1,u),v);
3327: case 3: return -kross(-1,u)*kross(-2*kross(-1,u),u*v);
3328: case 4: return v6==11 ? kross(-2,x1) : -kross(-2,u);
3329: }
3330: case -5:
3331: if (n2==1) return x1%32==23 ? 1 : -1;
3332: else return -kross(2,2*u+v);
3333: case -6:
3334: switch(n2)
3335: {
3336: case 1: return 1;
3337: case 2: return v6==10 ? 1 : -1;
3338: case 3: return (u%16==11) || ((u+4*v)%16==3) ? 1 : -1;
3339: }
3340: case -7:
3341: if (n2==1) return 1;
3342: else
3343: {
3344: y1=itos(modis(gaddsg(u,gmul2n(c6,-8)),16)); avma=av;
3345: if (v6==10) return (y1==9) || (y1==13) ? 1 : -1;
3346: else return (y1==9) || (y1==5) ? 1 : -1;
3347: }
3348: case -8: return n2==2 ? kross(-1,v*d1) : -1;
3349: case -9: return n2==2 ? -kross(-1,d1) : -1;
3350: default: return -1;
3351: }
3352: }
3353:
3354: static long
3355: ellrootno_3(GEN e)
3356: {
3357: long n2,kod,u,v,d1,av=avma,r6,k4,k6,v4;
3358: GEN p=stoi(3),c4,c6,tmp,p4;
3359:
3360: n2=neron(e,p,&kod); c4=(GEN)e[10]; c6=(GEN)e[11]; p4=stoi(81);
3361: if (gcmp0(c4)) { v4=12; u=0; }
3362: else { v4=pvaluation(c4,p,&tmp); u=itos(modii(tmp,p4)); }
3363: if (gcmp0(c6)) v=0;
3364: else {(void)pvaluation(c6,p,&tmp); v=itos(modii(tmp,p4));}
3365: (void)pvaluation((GEN)e[12],p,&tmp); d1=itos(modii(tmp,p4));
3366: avma=av;
3367: r6=v%9; k4=kross(u,3); k6=kross(v,3);
3368: if (kod>4) return k6;
3369: switch(kod)
3370: {
3371: case 1: case 3: case -3: return 1;
3372: case 2:
3373: switch(n2)
3374: {
3375: case 1: return (r6==4 || r6>6) ? 1 : -1;
3376: case 2: return -k4*k6;
3377: case 3: return 1;
3378: case 4: return -k6;
3379: }
3380: case 4:
3381: switch(n2)
3382: {
3383: case 1: return k6*kross(d1,3);
3384: case 2: return -k4;
3385: case 3: return -k6;
3386: }
3387: case -2: return n2==2 ? 1 : k6;
3388: case -4:
3389: switch(n2)
3390: {
3391: case 1:
3392: if (v4==4) return (r6==4 || r6==8) ? 1 : -1;
3393: else return (r6==1 || r6==2) ? 1 : -1;
3394: case 2: return -k6;
3395: case 3: return (r6==2 || r6==7) ? 1 : -1;
3396: case 4: return k6;
3397: }
3398: default: return -1;
3399: }
3400: }
3401:
3402: static long
3403: ellrootno_not23(GEN e, GEN p, GEN ex)
3404: {
3405: GEN j;
3406: long ep,z;
3407:
3408: if (gcmp1(ex)) return -kronecker(negi((GEN)e[11]),p);
3409: j=(GEN)e[13];
3410: if (!gcmp0(j) && ggval(j,p) < 0) return kronecker(negi(gun),p);
3411: ep=12/cgcd(12,ggval((GEN)e[12],p));
3412: if (ep==4) z=2;
3413: else z=(ep%2==0) ? 1 : 3;
3414: return kronecker(stoi(-z),p);
3415: }
3416:
3417: static long
3418: ellrootno_intern(GEN e, GEN p, GEN ex)
3419: {
3420: if (cmpis(p,3) > 0) return ellrootno_not23(e,p,ex);
3421: switch(itos(p))
3422: {
3423: case 3: return ellrootno_3(e);
3424: case 2: return ellrootno_2(e);
3425: default: err(talker,"incorrect prime in ellrootno_intern");
3426: }
3427: return 0; /* not reached */
3428: }
3429:
3430: /* local epsilon factor at p, including p=0 for the infinite place. Global
3431: if p==1. The equation can be non minimal, but must be over Q. Internal,
3432: no garbage collection. */
3433: static long
3434: ellrootno_all(GEN e, GEN p, GEN* ptcond)
3435: {
3436: long s,exs,i;
3437: GEN fa,gr,cond,pr,ex;
3438:
3439: gr=globalreduction(e);
3440: e=coordch(e,(GEN)gr[2]);
3441: cond=(GEN)gr[1]; if(ptcond) *ptcond=cond;
3442: if (typ(e[12]) != t_INT)
3443: err(talker,"not an integral curve in ellrootno");
3444: if (typ(p) != t_INT || signe(p)<0)
3445: err(talker,"not a nonnegative integer second arg in ellrootno");
3446: if (cmpis(p,2)>=0)
3447: {
3448: exs=ggval(cond,p);
3449: if (!exs) return 1;
3450: }
3451: if (cmpis(p,3)>0) return ellrootno_not23(e,p,stoi(exs));
3452: switch(itos(p))
3453: {
3454: case 3: return ellrootno_3(e);
3455: case 2: return ellrootno_2(e);
3456: case 1: s=-1; fa=factor(cond); pr=(GEN)fa[1]; ex=(GEN)fa[2];
3457: for (i=1; i<lg(pr); i++) s*=ellrootno_intern(e,(GEN)pr[i],(GEN)ex[i]);
3458: return s;
3459: case 0: return -1; /* local factor at infinity = -1 */
3460: default: return 0; /* never reached */
3461: }
3462: }
3463:
3464: long
3465: ellrootno(GEN e, GEN p)
3466: {
3467: long av=avma,s;
3468: if (!p) p = gun;
3469: s=ellrootno_all(e, p, NULL);
3470: avma=av; return s;
3471: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to elliptic.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / pari / src / modules