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