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