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