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