=================================================================== RCS file: /home/cvs/OpenXM_contrib/pari-2.2/src/modules/Attic/elliptic.c,v retrieving revision 1.1 retrieving revision 1.2 diff -u -p -r1.1 -r1.2 --- OpenXM_contrib/pari-2.2/src/modules/Attic/elliptic.c 2001/10/02 11:17:11 1.1 +++ OpenXM_contrib/pari-2.2/src/modules/Attic/elliptic.c 2002/09/11 07:27:05 1.2 @@ -1,4 +1,4 @@ -/* $Id: elliptic.c,v 1.1 2001/10/02 11:17:11 noro Exp $ +/* $Id: elliptic.c,v 1.2 2002/09/11 07:27:05 noro Exp $ Copyright (C) 2000 The PARI group. @@ -26,12 +26,11 @@ checkpt(GEN z) if (typ(z)!=t_VEC) err(elliper1); } -long +void checkell(GEN e) { long lx=lg(e); if (typ(e)!=t_VEC || lx<14) err(elliper1); - return lx; } void @@ -142,7 +141,7 @@ smallinitell0(GEN x, GEN y) GEN smallinitell(GEN x) { - ulong av = avma; + gpmem_t av = avma; GEN y = cgetg(14,t_VEC); smallinitell0(x,y); return gerepilecopy(av,y); } @@ -162,7 +161,7 @@ ellinit0(GEN x, long flag,long prec) void ellprint(GEN e) { - long av = avma; + gpmem_t av = avma; long vx = fetch_var(); long vy = fetch_var(); GEN z = cgetg(3,t_VEC); @@ -170,7 +169,7 @@ ellprint(GEN e) err(talker, "not an elliptic curve in ellprint"); z[1] = lpolx[vx]; name_var(vx, "X"); z[2] = lpolx[vy]; name_var(vy, "Y"); - fprintferr("%Z = %Z\n", ellLHS(e, z), ellRHS(e, polx[vx])); + fprintferr("%Z - (%Z)\n", ellLHS(e, z), ellRHS(e, polx[vx])); (void)delete_var(); (void)delete_var(); avma = av; } @@ -204,19 +203,21 @@ static GEN do_padic_agm(GEN *ptx1, GEN a1, GEN b1, GEN p) { GEN p1,r1,a,b,x,bmod1, bmod = modii((GEN)b1[4],p), x1 = *ptx1; - + long mi; + if (!x1) x1 = gmul2n(gsub(a1,b1),-2); + mi = min(precp(a1),precp(b1)); for(;;) { a=a1; b=b1; x=x1; - b1=gsqrt(gmul(a,b),0); bmod1=modii((GEN)b1[4],p); + b1=gprec(gsqrt(gmul(a,b),0),mi); bmod1=modii((GEN)b1[4],p); if (!egalii(bmod1,bmod)) b1 = gneg_i(b1); - a1=gmul2n(gadd(gadd(a,b),gmul2n(b1,1)),-2); + a1=gprec(gmul2n(gadd(gadd(a,b),gmul2n(b1,1)),-2),mi); r1=gsub(a1,b1); + if (gcmp0(r1)) {x1=x; break;} p1=gsqrt(gdiv(gadd(x,r1),x),0); if (! gcmp1(modii((GEN)p1[4],p))) p1 = gneg_i(p1); x1=gmul(x,gsqr(gmul2n(gaddsg(1,p1),-1))); - if (gcmp0(r1)) break; } *ptx1 = x1; return ginv(gmul2n(a1,2)); } @@ -227,13 +228,17 @@ padic_initell(GEN y, GEN p, long prec) GEN b2,b4,c4,c6,p1,p2,w,pv,a1,b1,x1,u2,q,e0,e1; long i,alpha; - if (valp(y[13]) >= 0) /* p | j */ + q=gadd(gun,ggrandocp(p,prec)); + for (i=1; i<=13; i++) y[i]=lmul(q,(GEN)y[i]); + if (gcmp0((GEN)y[13]) || valp((GEN)y[13]) >= 0) /* p | j */ err(talker,"valuation of j must be negative in p-adic ellinit"); if (egalii(p,gdeux)) - err(impl,"initell for 2-adic numbers"); /* pv=stoi(4); */ + { + pv=stoi(4); + err(impl,"initell for 2-adic numbers"); + } + else pv=p; - pv=p; q=ggrandocp(p,prec); - for (i=1; i<=5; i++) y[i]=ladd(q,(GEN)y[i]); b2= (GEN)y[6]; b4= (GEN)y[7]; c4= (GEN)y[10]; @@ -275,6 +280,16 @@ padic_initell(GEN y, GEN p, long prec) y[19]=zero; return y; } +/* mis pour debugger do_padic_agm. On peut enlever quand on veut */ + +GEN +dopad(GEN a, GEN b, GEN pv) +{ + GEN x1=NULL; + + return ginv(gmul2n(do_padic_agm(&x1,a,b,pv),2)); +} + static int invcmp(GEN x, GEN y) { return -gcmp(x,y); } @@ -355,21 +370,19 @@ initell0(GEN x, long prec) GEN initell(GEN x, long prec) { - ulong av = avma; + gpmem_t av = avma; return gerepilecopy(av, initell0(x,prec)); } GEN -coordch(GEN e, GEN ch) +_coordch(GEN e, GEN u, GEN r, GEN s, GEN t) { - GEN y,p1,p2,v,v2,v3,v4,v6,r,s,t,u; - long i,lx = checkell(e); - ulong av = avma; + GEN y, p1, p2, v, v2, v3, v4, v6; + long i, lx = lg(e); + gpmem_t av = avma; - checkch(ch); - u=(GEN)ch[1]; r=(GEN)ch[2]; s=(GEN)ch[3]; t=(GEN)ch[4]; - y=cgetg(lx,t_VEC); - v=ginv(u); v2=gsqr(v); v3=gmul(v,v2);v4=gsqr(v2); v6=gsqr(v3); + y = cgetg(lx,t_VEC); + v = ginv(u); v2 = gsqr(v); v3 = gmul(v,v2);v4 = gsqr(v2); v6 = gsqr(v3); y[1] = lmul(v,gadd((GEN)e[1],gmul2n(s,1))); y[2] = lmul(v2,gsub(gadd((GEN)e[2],gmulsg(3,r)),gmul(s,gadd((GEN)e[1],s)))); p2 = ellLHS0(e,r); @@ -388,42 +401,44 @@ coordch(GEN e, GEN ch) y[11] = lmul(v6,(GEN)e[11]); y[12] = lmul(gsqr(v6),(GEN)e[12]); y[13] = e[13]; - if (lx>14) + if (lx > 14) { - p1=(GEN)e[14]; + p1 = (GEN)e[14]; if (gcmp0(p1)) { y[14] = y[15] = y[16] = y[17] = y[18] = y[19] = zero; } + else if (typ(e[1])==t_PADIC) + { + y[14] = (long)_vec( gmul(v2, gsub((GEN)p1[1],r)) ); + y[15] = lmul((GEN)e[15], gsqr(u)); + y[16] = lmul((GEN)e[16], u); + y[17] = e[17]; + y[18] = e[18]; /* FIXME: how do q and w change ? */ + y[19] = zero; + } else { - if (typ(e[1])==t_PADIC) - { - p2=cgetg(2,t_VEC); p2[1]=lmul(v2,gsub((GEN)p1[1],r)); - y[14]=(long)p2; - y[15]=lmul(gsqr(u),(GEN)e[15]); - y[16]=lmul(u,(GEN)e[16]); -/* FIXME: how do q and w change ??? */ - y[17]=e[17]; - y[18]=e[18]; - y[19]=zero; - } - else - { - p2=cgetg(4,t_COL); - for (i=1; i<=3; i++) p2[i]=lmul(v2,gsub((GEN)p1[i],r)); - y[14]=(long)p2; - y[15]=lmul(u,(GEN)e[15]); - y[16]=lmul(u,(GEN)e[16]); - y[17]=ldiv((GEN)e[17],u); - y[18]=ldiv((GEN)e[18],u); - y[19]=lmul(gsqr(u),(GEN)e[19]); - } + p2 = cgetg(4,t_COL); + for (i=1; i<=3; i++) p2[i] = lmul(v2, gsub((GEN)p1[i],r)); + y[14] = (long)p2; + y[15] = lmul((GEN)e[15], u); + y[16] = lmul((GEN)e[16], u); + y[17] = ldiv((GEN)e[17], u); + y[18] = ldiv((GEN)e[18], u); + y[19] = lmul((GEN)e[19], gsqr(u)); } } return gerepilecopy(av,y); } +GEN +coordch(GEN e, GEN ch) +{ + checkch(ch); checkell(e); + return _coordch(e, (GEN)ch[1], (GEN)ch[2], (GEN)ch[3], (GEN)ch[4]); +} + static GEN pointch0(GEN x, GEN v2, GEN v3, GEN mor, GEN s, GEN t) { @@ -431,9 +446,9 @@ pointch0(GEN x, GEN v2, GEN v3, GEN mor, GEN s, GEN t) if (lg(x) < 3) return x; - z = cgetg(3,t_VEC); p1=gadd((GEN)x[1],mor); - z[1] = lmul(v2,p1); - z[2] = lmul(v3,gsub((GEN)x[2],gadd(gmul(s,p1),t))); + z = cgetg(3,t_VEC); p1 = gadd((GEN)x[1],mor); + z[1] = lmul(v2, p1); + z[2] = lmul(v3, gsub((GEN)x[2], gadd(gmul(s,p1),t))); return z; } @@ -441,18 +456,23 @@ GEN pointch(GEN x, GEN ch) { GEN y,v,v2,v3,mor,r,s,t,u; - long tx,lx=lg(x),i; - ulong av = avma; + long tx, i, lx = lg(x); + gpmem_t av = avma; checkpt(x); checkch(ch); if (lx < 2) return gcopy(x); - u=(GEN)ch[1]; r=(GEN)ch[2]; s=(GEN)ch[3]; t=(GEN)ch[4]; - tx=typ(x[1]); v=ginv(u); v2=gsqr(v); v3=gmul(v,v2); mor=gneg_i(r); + u = (GEN)ch[1]; + r = (GEN)ch[2]; + s = (GEN)ch[3]; + t = (GEN)ch[4]; + v = ginv(u); v2 = gsqr(v); v3 = gmul(v,v2); + mor = gneg_i(r); + tx = typ(x[1]); if (is_matvec_t(tx)) { - y=cgetg(lx,tx); + y = cgetg(lx,tx); for (i=1; i>1; ss=gadd(ss,gmul((GEN)z2[2],gpuigs(polx[0],ep))); - z2=gsub(z2,gmul((GEN)z2[2],gpuigs(z1,ep))); + ep=(-valp(z2))>>1; ss=gadd(ss,gmul((GEN)z2[2],gpowgs(polx[0],ep))); + z2=gsub(z2,gmul((GEN)z2[2],gpowgs(z1,ep))); } while (valp(z2)<=0); p2=gadd(p0,gmul(ss,p1)); p0=p1; p1=p2; @@ -646,7 +669,8 @@ GEN powell(GEN e, GEN z, GEN n) { GEN y; - long av=avma,i,j,tetpil,s; + long i, j, s; + gpmem_t av=avma, tetpil; ulong m; checksell(e); checkpt(z); @@ -679,7 +703,7 @@ mathell(GEN e, GEN x, long prec) { GEN y,p1,p2, *pdiag; long lx=lg(x),i,j,tx=typ(x); - ulong av = avma; + gpmem_t av = avma; if (!is_vec_t(tx)) err(elliper1); lx=lg(x); y=cgetg(lx,t_MAT); pdiag=(GEN*) new_chunk(lx); @@ -704,7 +728,8 @@ mathell(GEN e, GEN x, long prec) static GEN bilhells(GEN e, GEN z1, GEN z2, GEN h2, long prec) { - long lz1=lg(z1),tx,av=avma,tetpil,i; + long lz1=lg(z1), tx, i; + gpmem_t av=avma, tetpil; GEN y,p1,p2; if (lz1==1) return cgetg(1,typ(z1)); @@ -726,7 +751,8 @@ GEN bilhell(GEN e, GEN z1, GEN z2, long prec) { GEN p1,h2; - long av=avma,tetpil,tz1=typ(z1),tz2=typ(z2); + long tz1=typ(z1), tz2=typ(z2); + gpmem_t av=avma, tetpil; if (!is_matvec_t(tz1) || !is_matvec_t(tz2)) err(elliper1); if (lg(z1)==1) return cgetg(1,tz1); @@ -775,7 +801,8 @@ new_coords(GEN e, GEN x, GEN *pta, GEN *ptb, long prec GEN zell(GEN e, GEN z, long prec) { - long av=avma,ty,sw,fl; + long ty, sw, fl; + gpmem_t av=avma; GEN t,u,p1,p2,a,b,x1,u2,D = (GEN)e[12]; checkbell(e); @@ -857,15 +884,23 @@ zell(GEN e, GEN z, long prec) return gerepileupto(av,t); } +typedef struct { + GEN w1,w2; /* original basis for L = */ + GEN W1,W2,tau; /* new basis for L = = W2 <1,tau> */ + GEN a,b,c,d; /* tau in F = h/Sl2, tau = g.t, g=[a,b;c,d] in SL(2,Z) */ + GEN x,y; /* z/w2 defined mod <1,tau> --> z + x tau + y reduced mod <1,tau> */ +} SL2_red; + /* compute gamma in SL_2(Z) and t'=gamma(t) so that t' is in the usual fundamental domain. Internal function no check, no garbage. */ -static GEN -getgamma(GEN *ptt) +static void +set_gamma(SL2_red *T) { - GEN t = *ptt,a,b,c,d,n,m,p1,p2,run; + GEN t = T->tau,a,b,c,d,n,m,p1,run; - run = gsub(realun(DEFAULTPREC), gpuigs(stoi(10),-8)); - a=d=gun; b=c=gzero; + run = gsub(realun(DEFAULTPREC), gpowgs(stoi(10),-8)); + a = d = gun; + b = c = gzero; for(;;) { n = ground(greal(t)); @@ -880,142 +915,182 @@ getgamma(GEN *ptt) p1=negi(c); c=a; a=p1; p1=negi(d); d=b; b=p1; } - m=cgetg(3,t_MAT); *ptt = t; - p1=cgetg(3,t_COL); m[1]=(long)p1; - p2=cgetg(3,t_COL); m[2]=(long)p2; - p1[1]=(long)a; p2[1]=(long)b; - p1[2]=(long)c; p2[2]=(long)d; return m; + T->a = a; T->b = b; + T->c = c; T->d = d; } -static GEN -get_tau(GEN *ptom1, GEN *ptom2, GEN *ptga) +#define swap(x,y) { GEN _t=x; x=y; y=_t; } + +/* swap w1, w2 so that Im(t := w1/w2) > 0. Set tau = representative of t in + * the standard fundamental domain, and g in Sl_2, such that tau = g.t */ +static void +red_modSL2(SL2_red *T) { - GEN om1 = *ptom1, om2 = *ptom2, tau = gdiv(om1,om2); - long s = gsigne(gimag(tau)); - if (!s) - err(talker,"omega1 and omega2 R-linearly dependent in elliptic function"); - if (s < 0) { *ptom1=om2; *ptom2=om1; tau=ginv(tau); } - *ptga = getgamma(&tau); return tau; + long s; + T->tau = gdiv(T->w1,T->w2); + s = gsigne(gimag(T->tau)); + if (!s) err(talker,"w1 and w2 R-linearly dependent in elliptic function"); + if (s < 0) { swap(T->w1, T->w2); T->tau=ginv(T->tau); } + set_gamma(T); + /* update lattice */ + T->W1 = gadd(gmul(T->a,T->w1), gmul(T->b,T->w2)); + T->W2 = gadd(gmul(T->c,T->w1), gmul(T->d,T->w2)); + T->tau = gdiv(T->W1, T->W2); } static int -get_periods(GEN e, GEN *om1, GEN *om2) +get_periods(GEN e, SL2_red *T) { long tx = typ(e); if (is_vec_t(tx)) switch(lg(e)) { - case 3: *om1=(GEN)e[1]; *om2=(GEN)e[2]; return 1; - case 20: *om1=(GEN)e[16]; *om2=(GEN)e[15]; return 1; + case 3: T->w1=(GEN)e[1]; T->w2=(GEN)e[2]; red_modSL2(T); return 1; + case 20: T->w1=(GEN)e[16]; T->w2=(GEN)e[15];red_modSL2(T); return 1; } return 0; } -extern GEN PiI2(long prec); +static GEN +check_real(GEN q) +{ + return (typ(q) == t_COMPLEX && gcmp0((GEN)q[2]))? (GEN)q[1]: q; +} -/* computes the numerical values of eisenstein series. k is equal to a positive - even integer. If k=4 or 6, computes g2 or g3. If k=2, or k>6 even, - compute (2iPi/om2)^k*(1+2/zeta(1-k)*sum(n>=1,n^(k-1)q^n/(1-q^n)) with no - constant factor. */ GEN -elleisnum(GEN om, long k, long flag, long prec) +_elleisnum(SL2_red *T, long k, long prec) { - long av=avma,lim,av1; - GEN om1,om2,p1,pii2,tau,q,y,qn,ga,court,asub = NULL; /* gcc -Wall */ + gpmem_t lim, av; + GEN p1,pii2,q,y,qn,n; - if (k%2 || k<=0) err(talker,"k not a positive even integer in elleisnum"); - if (!get_periods(om, &om1, &om2)) err(typeer,"elleisnum"); pii2 = PiI2(prec); - tau = get_tau(&om1,&om2, &ga); - if (k==2) asub=gdiv(gmul(pii2,mulsi(12,gcoeff(ga,2,1))),om2); - om2=gadd(gmul(gcoeff(ga,2,1),om1),gmul(gcoeff(ga,2,2),om2)); - if (k==2) asub=gdiv(asub,om2); - q=gexp(gmul(pii2,tau),prec); - y=gzero; court=stoi(3); - av1=avma; lim=stack_lim(av1,1); qn=gun; court[2]=0; + q = gexp(gmul(pii2,T->tau), prec); + q = check_real(q); + y = gzero; n = utoi(1); + av = avma; lim = stack_lim(av,1); qn = gun; n[2] = 0; for(;;) { - court[2]++; qn=gmul(q,qn); - p1=gdiv(gmul(gpuigs(court,k-1),qn),gsub(gun,qn)); - y=gadd(y,p1); + n[2]++; qn = gmul(q,qn); + p1 = gdiv(gmul(gpowgs(n,k-1),qn), gsub(gun,qn)); if (gcmp0(p1) || gexpo(p1) <= - bit_accuracy(prec) - 5) break; - if (low_stack(lim, stack_lim(av1,1))) + y = gadd(y, p1); + if (low_stack(lim, stack_lim(av,1))) { GEN *gptr[2]; gptr[0]=&y; gptr[1]=&qn; if(DEBUGMEM>1) err(warnmem,"elleisnum"); - gerepilemany(av1,gptr,2); + gerepilemany(av,gptr,2); } } + y = gadd(gun, gmul(y, gdiv(gdeux, gzeta(stoi(1-k),prec)))); /* = E_k(tau) */ - y=gadd(gun,gmul(gdiv(gdeux,gzeta(stoi(1-k),prec)),y)); - p1=gpuigs(gdiv(pii2,om2),k); - y = gmul(p1,y); - if (k==2) y=gsub(y,asub); - else if (k==4 && flag) y=gdivgs(y,12); - else if (k==6 && flag) y=gdivgs(y,216); - return gerepileupto(av,y); + y = gmul(y, gpowgs(gdiv(pii2,T->W2),k)); + return check_real(y); } -/* compute eta1, eta2 */ +/* Return (2iPi)^k E_k(L) = (2iPi/w2)^k E_k(tau), with L = , k > 0 even + * E_k(tau) = 1 + 2/zeta(1-k) * sum(n>=1, n^(k-1) q^n/(1-q^n)) + * If flag is != 0 and k=4 or 6, compute g2 = E4/12 or g3 = E6/216 resp. */ GEN -elleta(GEN om, long prec) +elleisnum(GEN om, long k, long flag, long prec) { - long av=avma; + gpmem_t av = avma; + GEN p1, y; + SL2_red T; + + if (k&1 || k<=0) err(talker,"k not a positive even integer in elleisnum"); + if (!get_periods(om, &T)) err(typeer,"elleisnum"); + y = _elleisnum(&T, k, prec); + if (k==2 && signe(T.c)) + { + p1 = gmul(PiI2(prec), mulsi(12, T.c)); + y = gsub(y, gdiv(p1, gmul(T.w2, T.W2))); + } + else if (k==4 && flag) y = gdivgs(y, 12); + else if (k==6 && flag) y = gdivgs(y,216); + return gerepileupto(av,y); +} + +static GEN +_elleta(SL2_red *T, long prec) +{ GEN e2,y1,y2,y; - e2 = gdivgs(elleisnum(om,2,0,prec),12); - y2 = gmul((GEN)om[2],e2); - y1 = gadd(gdiv(PiI2(prec),(GEN)om[2]), gmul((GEN)om[1],e2)); + e2 = gdivgs(_elleisnum(T,2,prec),12); + y2 = gmul(T->W2, e2); + y1 = gadd(gdiv(PiI2(prec), T->W2), gmul(T->W1,e2)); y = cgetg(3,t_VEC); y[1] = lneg(y1); - y[2] = lneg(y2); return gerepileupto(av, y); + y[2] = lneg(y2); return y; } -/* computes the numerical value of wp(z | om1 Z + om2 Z), - If flall=1, compute also wp'. Reduce to the fundamental domain first. */ +/* compute eta1, eta2 */ +GEN +elleta(GEN om, long prec) +{ + gpmem_t av = avma; + SL2_red T; + if (!get_periods(om, &T)) err(typeer,"elleta"); + return gerepileupto(av, _elleta(&T,prec)); +} + static GEN -weipellnumall(GEN om1, GEN om2, GEN z, long flall, long prec) +reduce_z(GEN z, long prec, SL2_red *T) { - long av=avma,tetpil,lim,av1,toadd; - GEN p1,pii2,pii4,a,tau,q,u,y,yp,u1,u2,qn,v,ga; + GEN Z = gdiv(z, T->W2); + long t = typ(z); - pii2 = PiI2(prec); - tau = get_tau(&om1,&om2, &ga); - om2=gadd(gmul(gcoeff(ga,2,1),om1),gmul(gcoeff(ga,2,2),om2)); - z=gdiv(z,om2); - a=ground(gdiv(gimag(z),gimag(tau))); z=gsub(z,gmul(a,tau)); - a=ground(greal(z)); z=gsub(z,a); - if (gcmp0(z) || gexpo(z) < 5 - bit_accuracy(prec)) - { - avma=av; v=cgetg(2,t_VEC); v[1]=zero; return v; - } + if (!is_scalar_t(t) || t == t_INTMOD || t == t_PADIC || t == t_POLMOD) + err(typeer,"reduction mod SL2 (reduce_z)"); + T->x = ground(gdiv(gimag(Z), gimag(T->tau))); + Z = gsub(Z, gmul(T->x,T->tau)); + T->y = ground(greal(Z)); + Z = gsub(Z, T->y); + if (gcmp0(Z) || gexpo(Z) < 5 - bit_accuracy(prec)) Z = NULL; /* z in L */ + return Z; +} - q=gexp(gmul(pii2,tau),prec); - u=gexp(gmul(pii2,z),prec); - u1=gsub(gun,u); u2=gsqr(u1); - y=gadd(gdivgs(gun,12),gdiv(u,u2)); - if (flall) yp=gdiv(gadd(gun,u),gmul(u1,u2)); - toadd=(long)ceil(9.065*gtodouble(gimag(z))); +/* computes the numerical value of wp(z | L), L = om1 Z + om2 Z + * return NULL if z in L. If flall=1, compute also wp' */ +static GEN +weipellnumall(SL2_red *T, GEN z, long flall, long prec) +{ + long toadd; + gpmem_t av=avma, lim, av1; + GEN p1,pii2,q,u,y,yp,u1,u2,qn,v; - av1=avma; lim=stack_lim(av1,1); qn=q; + z = reduce_z(z,prec, T); + if (!z) return NULL; + + /* Now L,z normalized to <1,tau>. z in fund. domain of <1, tau> */ + pii2 = PiI2(prec); + q = gexp(gmul(pii2,T->tau),prec); + u = gexp(gmul(pii2,z),prec); + u1= gsub(gun,u); u2=gsqr(u1); + y = gadd(ginv(utoi(12)), gdiv(u,u2)); + if (flall) yp = gdiv(gadd(gun,u), gmul(u1,u2)); + toadd = (long)ceil(9.065*gtodouble(gimag(z))); + + av1 = avma; lim = stack_lim(av1,1); qn = q; for(;;) { - GEN p2,qnu,qnu1,qnu2,qnu3,qnu4; + GEN qnu,qnu1,qnu2,qnu3,qnu4; - qnu=gmul(qn,u); qnu1=gsub(gun,qnu); qnu2=gsqr(qnu1); - qnu3=gsub(qn,u); qnu4=gsqr(qnu3); - p1=gsub(gmul(u,gadd(ginv(qnu2),ginv(qnu4))), - gmul2n(ginv(gsqr(gsub(gun,qn))),1)); - p1=gmul(qn,p1); - y=gadd(y,p1); + qnu = gmul(qn,u); /* q^n u */ + qnu1 = gsub(gun,qnu); /* 1 - q^n u */ + qnu2 = gsqr(qnu1); /* (1 - q^n u)^2 */ + qnu3 = gsub(qn,u); /* q^n - u */ + qnu4 = gsqr(qnu3); /* (q^n - u)^2 */ + p1 = gsub(gmul(u, gadd(ginv(qnu2),ginv(qnu4))), + gmul2n(ginv(gsqr(gsub(gun,qn))), 1)); + y = gadd(y, gmul(qn,p1)); if (flall) { - p2=gadd(gdiv(gadd(gun,qnu),gmul(qnu1,qnu2)), - gdiv(gadd(qn,u),gmul(qnu3,qnu4))); - p2=gmul(qn,p2); - yp=gadd(yp,p2); + p1 = gadd(gdiv(gadd(gun,qnu),gmul(qnu1,qnu2)), + gdiv(gadd(qn,u),gmul(qnu3,qnu4))); + + yp = gadd(yp, gmul(qn,p1)); } - qn=gmul(q,qn); + qn = gmul(q,qn); if (gexpo(qn) <= - bit_accuracy(prec) - 5 - toadd) break; if (low_stack(lim, stack_lim(av1,1))) { @@ -1025,44 +1100,43 @@ weipellnumall(GEN om1, GEN om2, GEN z, long flall, lon } } - pii2=gdiv(pii2,om2); - pii4=gsqr(pii2); - y = gmul(pii4,y); - if (flall) yp=gmul(u,gmul(gmul(pii4,pii2),yp)); - tetpil=avma; - if (flall) { v=cgetg(3,t_VEC); v[1]=lcopy(y); v[2]=lmul2n(yp,-1); } - else v=gcopy(y); - return gerepile(av,tetpil,v); + u1 = gdiv(pii2, T->W2); + u2 = gsqr(u1); + y = gmul(u2,y); /* y *= (2i pi / w2)^2 */ + if (flall) + { + yp = gmul(u, gmul(gmul(u1,u2),yp));/* yp *= u (2i pi / w2)^3 */ + v = cgetg(3,t_VEC); + v[1] = (long)y; + v[2] = lmul2n(yp,-1); + } + else v = y; + return gerepilecopy(av, v); } GEN ellzeta(GEN om, GEN z, long prec) { - long av=avma,tetpil,lim,av1,toadd; - GEN zinit,om1,om2,p1,pii2,tau,q,u,y,u1,qn,ga,x1,x2,et; + long toadd; + gpmem_t av=avma, lim, av1; + GEN Z,p1,pii2,q,u,y,u1,qn,et; + SL2_red T; - if (!get_periods(om, &om1, &om2)) err(typeer,"ellzeta"); - pii2 = PiI2(prec); - tau = get_tau(&om1,&om2, &ga); - om2=gadd(gmul(gcoeff(ga,2,1),om1),gmul(gcoeff(ga,2,2),om2)); - om1=gmul(tau,om2); om=cgetg(3,t_VEC); om[1]=(long)om1; om[2]=(long)om2; - z=gdiv(z,om2); + if (!get_periods(om, &T)) err(typeer,"ellzeta"); + Z = reduce_z(z, prec, &T); - x1=ground(gdiv(gimag(z),gimag(tau))); z=gsub(z,gmul(x1,tau)); - x2=ground(greal(z)); z=gsub(z,x2); zinit=gmul(z,om2); - et=elleta(om,prec); - et=gadd(gmul(x1,(GEN)et[1]),gmul(x2,(GEN)et[2])); - if (gcmp0(z) || gexpo(z) < 5 - bit_accuracy(prec)) - { - p1=ginv(zinit); tetpil=avma; return gerepile(av,tetpil,gadd(p1,et)); - } - q=gexp(gmul(pii2,tau),prec); - u=gexp(gmul(pii2,z),prec); + et = _elleta(&T,prec); + et = gadd(gmul(T.x,(GEN)et[1]), gmul(T.y,(GEN)et[2])); + if (!Z) return gerepileupto(av, gadd(ginv(z),et)); /* FIXME ??? */ + + pii2 = PiI2(prec); + q=gexp(gmul(pii2,T.tau),prec); + u=gexp(gmul(pii2,Z),prec); u1=gsub(u,gun); - y=gdiv(gmul(gsqr(om2),elleisnum(om,2,0,prec)),pii2); - y=gadd(ghalf,gdivgs(gmul(z,y),-12)); + y=gdiv(gmul(gsqr(T.W2),_elleisnum(&T,2,prec)),pii2); + y=gadd(ghalf,gdivgs(gmul(Z,y),-12)); y=gadd(y,ginv(u1)); - toadd=(long)ceil(9.065*gtodouble(gimag(z))); + toadd=(long)ceil(9.065*gtodouble(gimag(Z))); av1=avma; lim=stack_lim(av1,1); qn=q; for(;;) { @@ -1079,53 +1153,50 @@ ellzeta(GEN om, GEN z, long prec) } } - y=gmul(gdiv(pii2,om2),y); - tetpil=avma; - return gerepile(av,tetpil,gadd(y,et)); + y=gmul(gdiv(pii2,T.W2),y); + return gerepileupto(av, gadd(y,et)); } /* if flag=0, return ellsigma, otherwise return log(ellsigma) */ GEN -ellsigma(GEN om, GEN z, long flag, long prec) +ellsigma(GEN w, GEN z, long flag, long prec) { - long av=avma,lim,av1,toadd; - GEN zinit,om1,om2,p1,pii2,tau,q,u,y,y1,u1,qn,ga,negu,uinv,x1,x2,et,etnew,uhalf; + long toadd; + gpmem_t av=avma, lim, av1; + GEN Z,zinit,p1,pii2,q,u,y,y1,u1,qn,negu,uinv,et,etnew,uhalf; int doprod = (flag >= 2); int dolog = (flag & 1); + SL2_red T; - if (!get_periods(om, &om1, &om2)) err(typeer,"ellsigmaprod"); - pii2 = PiI2(prec); - tau = get_tau(&om1,&om2, &ga); - om2=gadd(gmul(gcoeff(ga,2,1),om1),gmul(gcoeff(ga,2,2),om2)); - om1=gmul(tau,om2); om=cgetg(3,t_VEC); om[1]=(long)om1; om[2]=(long)om2; - z=gdiv(z,om2); + if (!get_periods(w, &T)) err(typeer,"ellsigma"); + Z = reduce_z(z, prec, &T); - x1=ground(gdiv(gimag(z),gimag(tau))); z=gsub(z,gmul(x1,tau)); - x2=ground(greal(z)); z=gsub(z,x2); zinit=gmul(z,om2); - et=elleta(om,prec); - etnew=gadd(gmul(x1,(GEN)et[1]),gmul(x2,(GEN)et[2])); - etnew=gmul(etnew,gadd(gmul2n(gadd(gmul(x1,om1),gmul(x2,om2)),-1),zinit)); - if (mpodd(x1) || mpodd(x2)) etnew=gadd(etnew,gmul2n(pii2,-1)); - if (gexpo(z) < 5 - bit_accuracy(prec)) - { - if (dolog) - return gerepileupto(av, gadd(etnew,glog(zinit,prec))); - else - return gerepileupto(av, gmul(gexp(etnew,prec),zinit)); + et = _elleta(&T, prec); + etnew = gadd(gmul(T.x,(GEN)et[1]), gmul(T.y,(GEN)et[2])); + zinit = Z? gmul(Z,T.W2): gzero; + p1 = gadd(zinit, gmul2n(gadd(gmul(T.x,T.W1),gmul(T.y,T.W2)),-1)); + etnew = gmul(etnew, p1); + pii2 = PiI2(prec); + if (mpodd(T.x) || mpodd(T.y)) etnew = gadd(etnew, gmul2n(pii2,-1)); + if (!Z) + { /* FIXME ??? */ + if (dolog) Z = gadd(etnew, glog(z,prec)); + else Z = gmul(gexp(etnew,prec), z); + return gerepileupto(av, Z); } - y1 = gadd(etnew,gmul2n(gmul(gmul(z,zinit),(GEN)et[2]),-1)); + y1 = gadd(etnew,gmul2n(gmul(gmul(Z,zinit),(GEN)et[2]),-1)); /* 9.065 = 2*Pi/log(2) */ - toadd = (long)ceil(9.065*fabs(gtodouble(gimag(z)))); - uhalf = gexp(gmul(gmul2n(pii2,-1),z),prec); + toadd = (long)ceil(9.065*fabs(gtodouble(gimag(Z)))); + uhalf = gexp(gmul(gmul2n(pii2,-1),Z),prec); u = gsqr(uhalf); if (doprod) { /* use product */ - q=gexp(gmul(pii2,tau),prec); + q=gexp(gmul(pii2,T.tau),prec); uinv=ginv(u); u1=gsub(uhalf,ginv(uhalf)); - y=gdiv(gmul(om2,u1),pii2); + y=gdiv(gmul(T.W2,u1),pii2); av1=avma; lim=stack_lim(av1,1); qn=q; negu=stoi(-1); for(;;) @@ -1147,8 +1218,8 @@ ellsigma(GEN om, GEN z, long flag, long prec) { /* use sum */ GEN q8,qn2,urn,urninv; long n; - q8=gexp(gmul2n(gmul(pii2,tau),-3),prec); - q=gpuigs(q8,8); + q8=gexp(gmul2n(gmul(pii2,T.tau),-3),prec); + q=gpowgs(q8,8); u=gneg_i(u); uinv=ginv(u); y=gzero; av1=avma; lim=stack_lim(av1,1); qn=q; qn2=gun; @@ -1169,7 +1240,7 @@ ellsigma(GEN om, GEN z, long flag, long prec) } } - p1=gmul(q8,gmul(gdiv(gdiv((GEN)om[2],pii2),gpuigs(trueeta(tau,prec),3)),y)); + p1=gmul(q8,gmul(gdiv(gdiv(T.W2,pii2),gpowgs(trueeta(T.tau,prec),3)),y)); } if (dolog) @@ -1181,83 +1252,118 @@ ellsigma(GEN om, GEN z, long flag, long prec) GEN pointell(GEN e, GEN z, long prec) { - long av=avma,tetpil; - GEN y,yp,v,p1; + gpmem_t av = avma; + GEN v; + SL2_red T; - checkbell(e); - p1=weipellnumall((GEN)e[16],(GEN)e[15],z,1,prec); - if (lg(p1)==2) { avma=av; v=cgetg(2,t_VEC); v[1]=zero; return v; } - y = gsub((GEN)p1[1], gdivgs((GEN)e[6],12)); - yp = gsub((GEN)p1[2], gmul2n(ellLHS0(e,y),-1)); - tetpil=avma; v=cgetg(3,t_VEC); v[1]=lcopy(y); v[2]=lcopy(yp); - return gerepile(av,tetpil,v); + checkbell(e); (void)get_periods(e, &T); + v = weipellnumall(&T,z,1,prec); + if (!v) { avma = av; return _vec(gzero); } + v[1] = lsub((GEN)v[1], gdivgs((GEN)e[6],12)); + v[2] = lsub((GEN)v[2], gmul2n(ellLHS0(e,(GEN)v[1]),-1)); + return gerepilecopy(av, v); } -GEN -weipell(GEN e, long prec) +static GEN +_weipell(GEN c4, GEN c6, long PREC) { - long av1,tetpil,precres,i,k,l; - GEN res,p1,s,t; + long i, k, l, precres = 2*PREC; + gpmem_t av1, tetpil; + GEN P,p1,s,t, res = cgetg(precres+2,t_SER); - checkell(e); precres = 2*prec+2; - res=cgetg(precres,t_SER); res[1] = evalsigne(1) | evalvalp(-2) | evalvarn(0); - if (!prec) { setsigne(res,0); return res; } - for (i=3; i 3) + { + checkell(e); + c4 = (GEN)e[10]; + c6 = (GEN)e[11]; + } + else + { + c4 = elleisnum(e, 4, 0, prec); + c6 = elleisnum(e, 6, 0, prec); c6 = gneg(c6); + } + return _weipell(c4,c6,PREC); +} + +/* assume x a t_POL */ +int +is_simple_var(GEN x) +{ + return (degpol(x) == 1 && gcmp0((GEN)x[2]) && gcmp1((GEN)x[3])); +} + +GEN +ellwp0(GEN w, GEN z, long flag, long prec, long PREC) +{ + GEN v; + gpmem_t av = avma; + SL2_red T; + + if (!z) return weipell0(w,prec,PREC); if (typ(z)==t_POL) { - if (lgef(z) != 4 || !gcmp0((GEN)z[2]) || !gcmp1((GEN)z[3])) - err(talker,"expecting a simple variable in ellwp"); - v = weipell(om,PREC); setvarn(v, varn(z)); + if (!is_simple_var(z)) err(talker,"expecting a simple variable in ellwp"); + v = weipell0(w,prec,PREC); setvarn(v, varn(z)); return v; } - if (!get_periods(om, &om1, &om2)) err(typeer,"ellwp"); + if (!get_periods(w, &T)) err(typeer,"ellwp"); switch(flag) { - case 0: v=weipellnumall(om1,om2,z,0,prec); - if (typ(v)==t_VEC && lg(v)==2) { avma=av; v=gpuigs(z,-2); } + case 0: v = weipellnumall(&T,z,0,prec); + if (!v) { avma = av; v = gpowgs(z,-2); } return v; - case 1: v=weipellnumall(om1,om2,z,1,prec); - if (typ(v)==t_VEC && lg(v)==2) + case 1: v = weipellnumall(&T,z,1,prec); + if (!v) { - GEN p1 = gmul2n(gpuigs(z,3),1); - long tetpil=avma; - v=cgetg(3,t_VEC); - v[1]=lpuigs(z,-2); - v[2]=lneg(p1); return gerepile(av,tetpil,v); + GEN p1 = gmul2n(gpowgs(z,3),1); + gpmem_t tetpil = avma; + v = cgetg(3,t_VEC); + v[1] = lpuigs(z,-2); + v[2] = lneg(p1); return gerepile(av,tetpil,v); } return v; - case 2: return pointell(om,z,prec); + case 2: return pointell(w,z,prec); default: err(flagerr,"ellwp"); return NULL; } } @@ -1266,7 +1372,7 @@ ellwp0(GEN om, GEN z, long flag, long prec, long PREC) static GEN _a_2(GEN e) { - long av = avma; + gpmem_t av = avma; GEN unmodp = gmodulss(1,8); ulong e6 = itos((GEN)gmul(unmodp,(GEN)e[6])[2]); ulong e8 = itos((GEN)gmul(unmodp,(GEN)e[8])[2]); @@ -1282,7 +1388,8 @@ apell2_intern(GEN e, ulong p) if (p == 2) return _a_2(e); else { - ulong av = avma, i; + ulong i; + gpmem_t av = avma; GEN unmodp = gmodulss(1,p); ulong e6 = itos((GEN)gmul(unmodp,(GEN)e[6])[2]); ulong e8 = itos((GEN)gmul(unmodp,(GEN)e[8])[2]); @@ -1337,14 +1444,15 @@ multi_invmod(GEN x, GEN p) static GEN addsell(GEN e, GEN z1, GEN z2, GEN p) { - GEN p1,p2,x,x1,x2,y,y1,y2; - long av = avma; + GEN z,p1,p2,x,x1,x2,y,y1,y2; + gpmem_t av; if (!z1) return z2; if (!z2) return z1; x1 = (GEN)z1[1]; y1 = (GEN)z1[2]; x2 = (GEN)z2[1]; y2 = (GEN)z2[2]; + z = cgetg(3, t_VEC); av = avma; p2 = subii(x2, x1); if (p2 == gzero) { @@ -1356,11 +1464,12 @@ addsell(GEN e, GEN z1, GEN z2, GEN p) else p1 = subii(y2,y1); p1 = mulii(p1, mpinvmod(p2, p)); p1 = resii(p1, p); - x = subii(sqri(p1), addii(x1,x2)); x = modii(x,p); + x = subii(sqri(p1), addii(x1,x2)); y = negi(addii(y1, mulii(p1,subii(x,x1)))); - avma = av; p1 = cgetg(3,t_VEC); - p1[1] = licopy(x); - p1[2] = lmodii(y, p); return p1; + x = modii(x,p); + y = modii(y,p); avma = av; + z[1] = licopy(x); + z[2] = licopy(y); return z; } /* z1 <-- z1 + z2 */ @@ -1387,6 +1496,15 @@ addsell_part2(GEN e, GEN z1, GEN z2, GEN p, GEN p2inv) } static GEN +negsell(GEN f) +{ + GEN g = cgetg(3, t_VEC); + g[1] = f[1]; + g[2] = lnegi((GEN)f[2]); + return g; +} + +static GEN powsell(GEN e, GEN z, GEN n, GEN p) { GEN y; @@ -1396,9 +1514,8 @@ powsell(GEN e, GEN z, GEN n, GEN p) if (!s || !z) return NULL; if (s < 0) { - n = negi(n); y = cgetg(3,t_VEC); - y[2] = lnegi((GEN)z[2]); - y[1] = z[1]; z = y; + z = negsell(z); + n = negi(n); } if (is_pm1(n)) return z; @@ -1417,6 +1534,25 @@ powsell(GEN e, GEN z, GEN n, GEN p) return addsell(e,y,z,p); } +/* assume H.f = 0, return exact order of f */ +static GEN +exact_order(GEN H, GEN f, GEN cp4, GEN p) +{ + GEN P, e, g, h = H, fa = decomp(H); + long i, j, l; + + P = (GEN)fa[1]; l = lg(P); + e = (GEN)fa[2]; + for (i=1; i= k */ static void _fix(GEN x, long k) @@ -1425,68 +1561,77 @@ _fix(GEN x, long k) if (lgefint(y) < k) { GEN p1 = cgeti(k); affii(y,p1); *x = (long)p1; } } -/* low word of integer x */ -#define _low(x) (__x=(GEN)x, __x[lgefint(__x)-1]) - -/* compute a_p using Shanks/Mestre + Montgomery's trick. Assume p > 20, say */ +/* compute a_p using Shanks/Mestre + Montgomery's trick. Assume p > 457 */ GEN apell1(GEN e, GEN p) { - long *tx, *ty, *ti, av = avma, av2,pfinal,i,j,j2,s,flc,flcc,x,nb; - GEN p1,p2,p3,h,mfh,f,fh,fg,pordmin,u,v,p1p,p2p,acon,bcon,c4,c6,cp4,pts; - GEN __x; - - if (DEBUGLEVEL) timer2(); - p1 = gmodulsg(1,p); - c4 = gdivgs(gmul(p1,(GEN)e[10]), -48); - c6 = gdivgs(gmul(p1,(GEN)e[11]), -864); - pordmin = gceil(gmul2n(gsqrt(p,DEFAULTPREC),2)); - p1p = addsi(1,p); p2p = shifti(p1p,1); - x=0; flcc=0; flc = kronecker((GEN)c6[2],p); - u=c6; acon=gzero; bcon=gun; h=p1p; + long *tx, *ty, *ti, pfinal, i, j, s, k, k1, x, nb; + gpmem_t av = avma, av2; + GEN p1,h,mfh,f,fh,fg,pordmin,u,v,p1p,p2p,A,B,c4,c6,cp4,pts; tx = ty = ti = NULL; /* gcc -Wall */ + + if (DEBUGLEVEL) (void)timer2(); + c4 = gmod(gdivgs((GEN)e[10], -48), p); + c6 = gmod(gdivgs((GEN)e[11], -864), p); + pordmin = gceil(gmul2n(gsqrt(p,DEFAULTPREC),2)); /* ceil( sqrt(4p) ) */ + p1p = addsi(1,p); + p2p = shifti(p1p,1); + x = 0; u = c6; k1 = 0; k = kronecker(c6, p); + A = gzero; B = gun; h = p1p; for(;;) { - while (flc==flcc || !flc) + while (k == k1 || !k) { x++; - u = gadd(c6, gmulsg(x, gaddgs(c4,x*x))); - flc = kronecker((GEN)u[2],p); + u = modii(addii(c6, mulsi(x, addii(c4, mulss(x,x)))), p); + k = kronecker(u, p); } - flcc = flc; + k1 = k; + f = cgetg(3,t_VEC); - f[1] = (long)lift_intern(gmulsg(x,u)); - f[2] = (long)lift_intern(gsqr(u)); - cp4 = lift_intern(gmul(c4, (GEN)f[2])); + f[1] = lmodii(mulsi(x,u), p); + f[2] = lmodii(sqri(u), p); + cp4 = modii(mulii(c4, (GEN)f[2]), p); fh = powsell(cp4,f,h,p); - s = itos(gceil(gsqrt(gdiv(pordmin,bcon),DEFAULTPREC))) >> 1; - nb = min(128, s >> 1); + s = itos(gceil(gsqrt(gdiv(pordmin,B),DEFAULTPREC))) >> 1; /* look for h s.t f^h = 0 */ - if (bcon == gun) + if (B == gun) { /* first time: initialize */ tx = newbloc(s+1); ty = newbloc(s+1); ti = newbloc(s+1); } - else f = powsell(cp4,f,bcon,p); /* F */ + else f = powsell(cp4,f,B,p); /* F = B.f */ *tx = evaltyp(t_VECSMALL) | evallg(s+1); if (!fh) goto FOUND; p1 = gcopy(fh); - pts = new_chunk(nb+1); + if (s < 3) + { /* we're nearly done: naive search */ + GEN q1 = p1, mf = negsell(f); /* -F */ + for (i=1;; i++) + { + p1 = addsell(cp4,p1, f,p); /* h.f + i.F */ + if (!p1) { h = addii(h, mulsi( i,B)); goto FOUND; } + q1 = addsell(cp4,q1,mf,p); /* h.f - i.F */ + if (!q1) { h = addii(h, mulsi(-i,B)); goto FOUND; } + } + } + /* Baby Step/Giant Step */ + nb = min(128, s >> 1); /* > 0. Will do nb pts at a time: faster inverse */ + pts = cgetg(nb+1, t_VEC); j = lgefint(p); for (i=1; i<=nb; i++) { /* baby steps */ pts[i] = (long)p1; /* h.f + (i-1).F */ - _fix(p1+1, j); tx[i] = _low((GEN)p1[1]); - _fix(p1+2, j); ty[i] = _low((GEN)p1[2]); + _fix(p1+1, j); tx[i] = modBIL((GEN)p1[1]); + _fix(p1+2, j); ty[i] = modBIL((GEN)p1[2]); p1 = addsell(cp4,p1,f,p); /* h.f + i.F */ - if (!p1) { h = addii(h, mulsi(i,bcon)); goto FOUND; } + if (!p1) { h = addii(h, mulsi(i,B)); goto FOUND; } } - mfh = dummycopy(fh); - mfh[2] = lnegi((GEN)mfh[2]); - fg = addsell(cp4,p1,mfh,p); /* nb.F */ - if (!fg) { h = mulsi(nb,bcon); goto FOUND; } + mfh = negsell(fh); + fg = addsell(cp4,p1,mfh,p); /* h.f + nb.F - h.f = nb.F */ + if (!fg) { h = mulsi(nb,B); goto FOUND; } u = cgetg(nb+1, t_VEC); av2 = avma; /* more baby steps, nb points at a time */ while (i <= s) @@ -1499,9 +1644,8 @@ apell1(GEN e, GEN p) if (u[j] == zero) /* sum = 0 or doubling */ { long k = i+j-2; - if (egalii((GEN)p1[2],(GEN)fg[2])) k -= 2*nb; /* fg = p1 */ - h = addii(h, mulsi(k,bcon)); - goto FOUND; + if (egalii((GEN)p1[2],(GEN)fg[2])) k -= 2*nb; /* fg == p1 */ + h = addii(h, mulsi(k,B)); goto FOUND; } } v = multi_invmod(u, p); @@ -1510,18 +1654,19 @@ apell1(GEN e, GEN p) { p1 = (GEN)pts[j]; addsell_part2(cp4,p1,fg,p, (GEN)v[j]); - tx[i] = _low((GEN)p1[1]); - ty[i] = _low((GEN)p1[2]); + tx[i] = modBIL((GEN)p1[1]); + ty[i] = modBIL((GEN)p1[2]); } avma = av2; } - p1 = addsell(cp4,(GEN)pts[j-1],mfh,p); /* = f^(s-1) */ + p1 = addsell(cp4,(GEN)pts[j-1],mfh,p); /* = (s-1).F */ + if (!p1) { h = mulsi(s-1,B); goto FOUND; } if (DEBUGLEVEL) msgtimer("[apell1] baby steps, s = %ld",s); - /* giant steps: fg = f^s */ + /* giant steps: fg = s.F */ fg = addsell(cp4,p1,f,p); - if (!fg) { h = addii(h, mulsi(s,bcon)); goto FOUND; } - pfinal = _low(p); av2 = avma; + if (!fg) { h = mulsi(s,B); goto FOUND; } + pfinal = modBIL(p); av2 = avma; p1 = gen_sort(tx, cmp_IND | cmp_C, NULL); for (i=1; i<=s; i++) ti[i] = tx[p1[i]]; @@ -1531,8 +1676,12 @@ apell1(GEN e, GEN p) avma = av2; gaffect(fg, (GEN)pts[1]); - for (j=2; j<=nb; j++) /* pts = first nb multiples of fg */ - gaffect(addsell(cp4,(GEN)pts[j-1],fg,p), (GEN)pts[j]); + for (j=2; j<=nb; j++) /* pts[j] = j.fg = (s*j).F */ + { + p1 = addsell(cp4,(GEN)pts[j-1],fg,p); + if (!p1) { h = mulii(mulss(s,j), B); goto FOUND; } + gaffect(p1, (GEN)pts[j]); + } /* replace fg by nb.fg since we do nb points at a time */ avma = av2; fg = gcopy((GEN)pts[nb]); @@ -1542,10 +1691,10 @@ apell1(GEN e, GEN p) { GEN ftest = (GEN)pts[j]; ulong m, l = 1, r = s+1; - long k, k2; + long k, k2, j2; avma = av2; - k = _low((GEN)ftest[1]); + k = modBIL((GEN)ftest[1]); while (l> 1; @@ -1554,24 +1703,24 @@ apell1(GEN e, GEN p) if (r <= (ulong)s && tx[r] == k) { while (tx[r] == k && r) r--; - k2 = _low((GEN)ftest[2]); + k2 = modBIL((GEN)ftest[2]); for (r++; tx[r] == k && r <= (ulong)s; r++) if (ty[r] == k2 || ty[r] == pfinal - k2) { /* [h+j2] f == ± ftest (= [i.s] f)? */ - if (DEBUGLEVEL) msgtimer("[apell1] giant steps, i = %ld",i); j2 = ti[r] - 1; + if (DEBUGLEVEL) msgtimer("[apell1] giant steps, i = %ld",i); p1 = addsell(cp4, powsell(cp4,f,stoi(j2),p),fh,p); if (egalii((GEN)p1[1], (GEN)ftest[1])) { if (egalii((GEN)p1[2], (GEN)ftest[2])) i = -i; - h = addii(h, mulii(addis(mulss(s,i), j2), bcon)); + h = addii(h, mulii(addis(mulss(s,i), j2), B)); goto FOUND; } } } if (++j > nb) { /* compute next nb points */ - long save = 0; /* gcc -Wall */ + long save = 0; /* gcc -Wall */; for (j=1; j<=nb; j++) { p1 = (GEN)pts[j]; @@ -1590,35 +1739,29 @@ apell1(GEN e, GEN p) } } -FOUND: /* success, found a point of exponent h */ - p2 = decomp(h); p1=(GEN)p2[1]; p2=(GEN)p2[2]; - for (i=1; i bcon) acon -= bcon; + A = (p2p - A) % B; + if ((A << 1) > B) A -= B; } - q = ((ulong)(p2p + bcon - (acon << 1))) / (bcon << 1); - h = acon + q*bcon; + q = ((ulong)(p2p + B - (A << 1))) / (B << 1); + h = A + q*B; avma = av; if (!i) break; } - free(table); return stoi((flc==1)? p1p-h: h-p1p); + free(table); return stoi(k==1? p1p-h: h-p1p); } GEN @@ -1785,7 +1946,8 @@ apell(GEN e, GEN p) if (typ(p)!=t_INT || signe(p)<0) err(talker,"not a prime in apell"); if (gdivise((GEN)e[12],p)) /* e[12] may be an intmod */ { - long av = avma,s; + long s; + gpmem_t av = avma; GEN p0 = egalii(p,gdeux)? stoi(8): p; GEN c6 = gmul((GEN)e[11],gmodulsg(1,p0)); s = kronecker(lift_intern(c6),p); avma=av; @@ -1809,7 +1971,8 @@ GEN anell(GEN e, long n) { long tab[4]={0,1,1,-1}; /* p prime; (-1/p) = tab[p&3]. tab[0] is not used */ - long p, i, m, av, tetpil; + long p, i, m; + gpmem_t av, tetpil; GEN p1,p2,an; checkell(e); @@ -1873,7 +2036,8 @@ anell(GEN e, long n) GEN akell(GEN e, GEN n) { - long i,j,ex,av=avma; + long i, j, ex; + gpmem_t av=avma; GEN p1,p2,ap,u,v,w,y,pl; checkell(e); @@ -1908,7 +2072,8 @@ akell(GEN e, GEN n) GEN hell(GEN e, GEN a, long prec) { - long av=avma,tetpil,n; + long n; + gpmem_t av=avma, tetpil; GEN p1,p2,y,z,q,pi2surw,pi2isurw,qn,ps; checkbell(e); @@ -1951,16 +2116,25 @@ hells(GEN e, GEN x, long prec) return mu; } +/* [1,0,0,0] */ +static GEN +init_ch() { + GEN v = cgetg(5,t_VEC); + v[1] = un; v[2] = v[3] = v[4] = zero; return v; +} + GEN hell2(GEN e, GEN x, long prec) { GEN ep,e3,ro,p1,p2,mu,d,xp; - long av=avma,tetpil,lx,lc,i,j,tx; + long lx, lc, i, j, tx; + gpmem_t av=avma, tetpil; if (!oncurve(e,x)) err(heller1); d=(GEN)e[12]; ro=(GEN)e[14]; e3=(gsigne(d) < 0)?(GEN)ro[1]:(GEN)ro[3]; - p1=cgetg(5,t_VEC); p1[1]=un; p1[2]=laddgs(gfloor(e3),-1); - p1[3]=p1[4]=zero; ep=coordch(e,p1); xp=pointch(x,p1); + p1 = init_ch(); p1[2] = laddgs(gfloor(e3),-1); + ep = coordch(e,p1); + xp = pointch(x,p1); tx=typ(x[1]); lx=lg(x); if (!is_matvec_t(tx)) { @@ -2016,7 +2190,8 @@ hell0(GEN e, GEN z, long prec) static GEN ghell0(GEN e, GEN a, long flag, long prec) { - long av=avma,lx,i,n,n2,grandn,tx=typ(a); + long lx, i, n, n2, grandn, tx=typ(a); + gpmem_t av=avma; GEN p,p1,p2,x,y,z,phi2,psi2,psi3,logdep; checkbell(e); if (!is_matvec_t(tx)) err(elliper1); @@ -2106,7 +2281,8 @@ static long ellrootno_all(GEN e, GEN p, GEN* ptcond); GEN lseriesell(GEN e, GEN s, GEN A, long prec) { - long av=avma,av1,tetpil,lim,l,n,eps,flun; + long l, n, eps, flun; + gpmem_t av=avma, av1, tetpil, lim; GEN z,p1,p2,cg,cg1,v,cga,cgb,s2,ns,gs,N; if (!A) A = gun; @@ -2118,14 +2294,15 @@ lseriesell(GEN e, GEN s, GEN A, long prec) } flun = gcmp1(A) && gcmp1(s); eps = ellrootno_all(e,gun,&N); - if (flun && eps<0) { z=cgetr(prec); affsr(0,z); return z; } + if (flun && eps < 0) return realzero(prec); + cg1=mppi(prec); setexpo(cg1,2); cg=divrr(cg1,gsqrt(N,prec)); cga=gmul(cg,A); cgb=gdiv(cg,A); l=(long)((pariC2*(prec-2) + fabs(gtodouble(s)-1.)*log(rtodbl(cga))) / rtodbl(cgb)+1); v = anell(e, min((ulong)l,TEMPMAX)); s2 = ns = NULL; /* gcc -Wall */ - if (!flun) { s2=gsubsg(2,s); ns=gpui(cg,gsubgs(gmul2n(s,1),2),prec); } + if (!flun) { s2=gsubsg(2,s); ns=gpow(cg,gsubgs(gmul2n(s,1),2),prec); } z=gzero; if (typ(s)==t_INT) { @@ -2136,9 +2313,9 @@ lseriesell(GEN e, GEN s, GEN A, long prec) av1=avma; lim=stack_lim(av1,1); for (n=1; n<=l; n++) { - p1=gdiv(incgam4(s,gmulsg(n,cga),gs,prec),gpui(stoi(n),s,prec)); + p1=gdiv(incgam4(s,gmulsg(n,cga),gs,prec),gpow(stoi(n),s,prec)); p2=flun? p1: gdiv(gmul(ns,incgam(s2,gmulsg(n,cgb),prec)), - gpui(stoi(n),s2,prec)); + gpow(stoi(n),s2,prec)); if (eps<0) p2=gneg_i(p2); z = gadd(z, gmul(gadd(p1,p2), ((ulong)n<=TEMPMAX)? (GEN)v[n]: akell(e,stoi(n)))); @@ -2177,60 +2354,70 @@ lseriesell(GEN e, GEN s, GEN A, long prec) Uses Tate's algorithm (Anvers IV). Given the remarks at the bottom of page 46, the "long" algorithm is used for p < 4 only. */ static void cumule(GEN *vtotal, GEN *e, GEN u, GEN r, GEN s, GEN t); -static void cumule1(GEN *vtotal, GEN *e, GEN v2); +static void cumulev(GEN *vtotal, GEN u, GEN r, GEN s, GEN t); static GEN -localreduction_result(long av, long f, long kod, long c, GEN v) +localred_result(long f, long kod, long c, GEN v) { - long tetpil = avma; - GEN result = cgetg(5, t_VEC); - result[1] = lstoi(f); result[2] = lstoi(kod); - result[3] = lcopy(v); result[4] = lstoi(c); - return gerepile(av,tetpil, result); + GEN z = cgetg(5, t_VEC); + z[1] = lstoi(f); + z[2] = lstoi(kod); + z[3] = lcopy(v); + z[4] = lstoi(c); return z; } -/* ici, p != 2 et p != 3 */ +/* Here p > 3. e assumed integral */ static GEN -localreduction_carac_not23(GEN e, GEN p) +localred_carac_p(GEN e, GEN p, int minim) { - long av = avma, k, f, kod, c, nuj, nudelta; - GEN pk, p2k, a2prime, a3prime; - GEN p2, r = gzero, s = gzero, t = gzero, v; - GEN c4, c6, delta, unmodp, xun, tri, var, p4k, p6k; + long k, f, kod, c, nuj, nudelta; + GEN p2, v = init_ch(); + GEN c4, c6, delta, tri; - nudelta = ggval((GEN)e[12], p); - v = cgetg(5,t_VEC); v[1] = un; v[2] = v[3] = v[4] = zero; - nuj = gcmp0((GEN)e[13])? 0: - ggval((GEN)e[13], p); + c4 = (GEN)e[10]; + c6 = (GEN)e[11]; + delta = (GEN)e[12]; + nuj = gcmp0((GEN)e[13])? 0: - ggval((GEN)e[13], p); + nudelta = ggval(delta, p); k = (nuj > 0 ? nudelta - nuj : nudelta) / 12; - c4 = (GEN)e[10]; c6 = (GEN)e[11]; delta = (GEN)e[12]; - if (k > 0) /* modele non minimal */ + if (k <= 0) { - pk = gpuigs(p, k); - if (mpodd((GEN)e[1])) - s = shifti(subii(pk, (GEN)e[1]), -1); - else - s = negi(shifti((GEN)e[1], -1)); - p2k = sqri(pk); - p4k = sqri(p2k); - p6k = mulii(p4k, p2k); + if (minim) return v; + } + else + { /* model not minimal */ + GEN pk = gpowgs(p,k), p2k = sqri(pk), p4k = sqri(p2k), p6k = mulii(p4k,p2k); + GEN r, s, t; - a2prime = subii((GEN)e[2], mulii(s, addii((GEN)e[1], s))); - switch(smodis(a2prime, 3)) + s = negi((GEN)e[1]); + if (mpodd(s)) s = addii(s, pk); + s = shifti(s, -1); + + r = subii(mulii(s, addii((GEN)e[1], s)), (GEN)e[2]); /* - a_2' */ + switch(smodis(r, 3)) { - case 0: r = negi(divis(a2prime, 3)); break; - case 1: r = divis(subii(p2k, a2prime), 3); break; - case 2: r = negi(divis(addii(a2prime, p2k), 3)); break; + default: break; /* 0 */ + case 1: r = addii(r, p2k); break; + case 2: r = subii(r, p2k); break; } - a3prime = ellLHS0_i(e,r); - if (mpodd(a3prime)) - t = shifti(subii(mulii(pk, p2k), a3prime), -1); - else - t = negi(shifti(a3prime, -1)); - v[1] = (long)pk; v[2] = (long)r; v[3] = (long)s; v[4] = (long)t; + r = divis(r, 3); + + t = negi(ellLHS0_i(e,r)); /* - a_3' */ + if (mpodd(t)) t = addii(t, mulii(pk, p2k)); + t = shifti(t, -1); + + v[1] = (long)pk; + v[2] = (long)r; + v[3] = (long)s; + v[4] = (long)t; + if (minim) return v; + nudelta -= 12 * k; - c4 = divii(c4, p4k); c6 = divii(c6, p6k); - delta = divii(delta, sqri(p6k)); + c4 = diviiexact(c4, p4k); + c6 = diviiexact(c6, p6k); + delta = diviiexact(delta, sqri(p6k)); } + if (nuj > 0) switch(nudelta - nuj) { case 0: f = 1; kod = 4+nuj; /* Inu */ @@ -2244,64 +2431,59 @@ localreduction_carac_not23(GEN e, GEN p) break; case 6: f = 2; kod = -4-nuj; /* Inu* */ if (nuj & 1) - c = 3 + kronecker(divii(mulii(c6, delta),gpuigs(p, 9+nuj)), p); + c = 3 + kronecker(divii(mulii(c6, delta),gpowgs(p, 9+nuj)), p); else - c = 3 + kronecker(divii(delta, gpuigs(p, 6+nuj)), p); + c = 3 + kronecker(divii(delta, gpowgs(p, 6+nuj)), p); break; default: err(bugparier,"localred (nu_delta - nu_j != 0,6)"); return NULL; /* not reached */ } else switch(nudelta) { - case 0: f = 0; kod = 1; c = 1; break; /* I0, regulier */ + case 0: f = 0; kod = 1; c = 1; break; /* I0, regular */ case 2: f = 2; kod = 2; c = 1; break; /* II */ case 3: f = 2; kod = 3; c = 2; break; /* III */ case 4: f = 2; kod = 4; /* IV */ - c = 2 + kronecker(gdiv(mulis(c6, -6), sqri(p)), p); + c = 2 + kronecker(divii(mulis(c6, -6), sqri(p)), p); break; case 6: f = 2; kod = -1; /* I0* */ p2 = sqri(p); - unmodp = gmodulsg(1,p); - var = gmul(unmodp,polx[0]); - tri = gsub(gsqr(var),gmul(divii(gmulsg(3, c4), p2),unmodp)); - tri = gsub(gmul(tri, var), - gmul(divii(gmul2n(c6,1), mulii(p2,p)),unmodp)); - xun = gmodulcp(var,tri); - c = lgef(ggcd((GEN)(gsub(gpui(xun,p,0),xun))[2], tri)) - 2; + /* x^3 - 3c4/p^2 x - 2c6/p^3 */ + tri = coefs_to_pol(4, gun, gzero, + negi(divii(gmul2n(c6,1), mulii(p2,p))), + negi(divii(gmulsg(3, c4), p2))); + c = 1 + FpX_nbroots(tri, p); break; case 8: f = 2; kod = -4; /* IV* */ - c = 2 + kronecker(gdiv(mulis(c6,-6), gpuigs(p,4)), p); + c = 2 + kronecker(gdiv(mulis(c6,-6), gpowgs(p,4)), p); break; case 9: f = 2; kod = -3; c = 2; break; /* III* */ case 10: f = 2; kod = -2; c = 1; break; /* II* */ default: err(bugparier,"localred"); return NULL; /* not reached */ } - return localreduction_result(av,f,kod,c,v); + return localred_result(f, kod, c, v); } -/* renvoie a_{ k,l } avec les notations de Tate */ +/* return a_{ k,l } in Tate's notation */ static int aux(GEN ak, int p, int l) { - long av = avma, pl = p, res; - while (--l) pl *= p; - res = smodis(divis(ak, pl), p); + gpmem_t av = avma; + long res = smodis(divis(ak, u_pow(p, l)), p); avma = av; return res; } static int aux2(GEN ak, int p, GEN pl) { - long av = avma, res; - res = smodis(divii(ak, pl), p); - avma = av; - return res; + gpmem_t av = avma; + long res = smodis(divii(ak, pl), p); + avma = av; return res; } -/* renvoie le nombre de racines distinctes du polynome XXX + aXX + bX + c - * modulo p s'il y a une racine multiple, elle est renvoyee dans *mult - */ +/* number of distinct roots of X^3 + aX^2 + bX + c modulo p + * if there's a multiple root, put it un *mult */ static int numroots3(int a, int b, int c, int p, int *mult) { @@ -2317,7 +2499,7 @@ numroots3(int a, int b, int c, int p, int *mult) } } -/* idem pour aXX +bX + c */ +/* same for aX^2 +bX + c */ static int numroots2(int a, int b, int c, int p, int *mult) { @@ -2325,39 +2507,40 @@ numroots2(int a, int b, int c, int p, int *mult) else { *mult = a * b; return (b * b - a * c) % 3 ? 2 : 1; } } -/* ici, p1 = 2 ou p1 = 3 */ +/* p = 2 or 3 */ static GEN -localreduction_carac_23(GEN e, GEN p1) +localred_carac_23(GEN e, long p, int minim) { - long av = avma, p, c, nu, nudelta; - int a21, a42, a63, a32, a64, theroot, al, be, ga; - GEN pk, p2k, pk1, p4, p6; - GEN p2, p3, r = gzero, s = gzero, t = gzero, v; + long c, nu, nudelta; + int a21, a42, a63, a32, a64, theroot, al, be, ga, p2, p3, p4; + GEN pk, p2k, pk1; + GEN r, s, t, v; - nudelta = ggval((GEN)e[12], p1); - v = cgetg(5,t_VEC); v[1] = un; v[2] = v[3] = v[4] = zero; + nudelta = ggval((GEN)e[12], stoi(p)); + v = init_ch(); - for(;;) + for (;;) { if (!nudelta) - return localreduction_result(av, 0, 1, 1, v); - /* I0 */ - p = itos(p1); - if (!divise((GEN)e[6], p1)) + return localred_result(0, 1, 1, v); + /* I0 */ + if (smodis((GEN)e[6], p)) /* p \nmid b2 */ { if (smodis(negi((GEN)e[11]), p == 2 ? 8 : 3) == 1) - c = nudelta; + c = nudelta; else - c = 2 - (nudelta & 1); - return localreduction_result(av, 1, 4 + nudelta, c, v); + c = 2 - (nudelta & 1); + return localred_result(1, 4 + nudelta, c, v); } - /* Inu */ + /* Inu */ if (p == 2) { r = modis((GEN)e[4], 2); s = modis(addii(r, (GEN)e[2]), 2); - if (signe(r)) t = modis(addii(addii((GEN)e[4], (GEN)e[5]), s), 2); - else t = modis((GEN)e[5], 2); + if (signe(r)) + t = modis(addii(addii((GEN)e[4], (GEN)e[5]), s), 2); + else + t = modis((GEN)e[5], 2); } else /* p == 3 */ { @@ -2365,259 +2548,296 @@ localreduction_carac_23(GEN e, GEN p1) s = modis((GEN)e[1], 3); t = modis(ellLHS0_i(e,r), 3); } - cumule(&v, &e, gun, r, s, t); /* p | a1, a2, a3, a4 et a6 */ - p2 = stoi(p*p); - if (!divise((GEN)e[5], p2)) - return localreduction_result(av, nudelta, 2, 1, v); - /* II */ - p3 = stoi(p*p*p); - if (!divise((GEN)e[9], p3)) - return localreduction_result(av, nudelta - 1, 3, 2, v); - /* III */ - if (!divise((GEN)e[8], p3)) + cumule(&v, &e, gun, r, s, t); /* p | (a1, a2, a3, a4, a6) */ + p2 = p * p; + if (smodis((GEN)e[5], p2)) + return localred_result(nudelta, 2, 1, v); + /* II */ + p3 = p2 * p; + if (smodis((GEN)e[9], p3)) + return localred_result(nudelta - 1, 3, 2, v); + /* III */ + if (smodis((GEN)e[8], p3)) { - if (smodis((GEN)e[8], (p==2)? 32: 27) == p*p) - c = 3; + if (smodis((GEN)e[8], (p==2)? 32: 27) == p2) + c = 3; else - c = 1; - return localreduction_result(av, nudelta - 2, 4, c, v); + c = 1; + return localred_result(nudelta - 2, 4, c, v); } - /* IV */ + /* IV */ - /* now for the last five cases... */ - - if (!divise((GEN)e[5], p3)) + if (smodis((GEN)e[5], p3)) cumule(&v, &e, gun, gzero, gzero, p == 2? gdeux: modis((GEN)e[3], 9)); - /* p | a1, a2; p^2 | a3, a4; p^3 | a6 */ - a21 = aux((GEN)e[2], p, 1); a42 = aux((GEN)e[4], p, 2); + /* p | a1, a2; p^2 | a3, a4; p^3 | a6 */ + a21 = aux((GEN)e[2], p, 1); + a42 = aux((GEN)e[4], p, 2); a63 = aux((GEN)e[5], p, 3); switch (numroots3(a21, a42, a63, p, &theroot)) { case 3: - if (p == 2) - c = 1 + (a63 == 0) + ((a21 + a42 + a63) & 1); - else - c = 1 + (a63 == 0) + (((1 + a21 + a42 + a63) % 3) == 0) - + (((1 - a21 + a42 - a63) % 3) == 0); - return localreduction_result(av, nudelta - 4, -1, c, v); - /* I0* */ - case 2: /* calcul de nu */ - if (theroot) cumule(&v, &e, gun, stoi(theroot * p), gzero, gzero); - /* p | a1; p^2 | a2, a3; p^3 | a4; p^4 | a6 */ - nu = 1; - pk = p2; - p2k = stoi(p * p * p * p); - for(;;) - { - if (numroots2(al = 1, be = aux2((GEN)e[3], p, pk), - ga = -aux2((GEN)e[5], p, p2k), p, &theroot) == 2) - break; - if (theroot) cumule(&v, &e, gun, gzero, gzero, mulsi(theroot,pk)); - pk1 = pk; pk = mulsi(p, pk); p2k = mulsi(p, p2k); - nu++; - if (numroots2(al = a21, be = aux2((GEN)e[4], p, pk), - ga = aux2((GEN)e[5], p, p2k), p, &theroot) == 2) - break; - if (theroot) cumule(&v, &e, gun, mulsi(theroot, pk1), gzero, gzero); - p2k = mulsi(p, p2k); - nu++; - } - if (p == 2) - c = 4 - 2 * (ga & 1); - else - c = 3 + kross(be * be - al * ga, 3); - return localreduction_result(av, nudelta - 4 - nu, -4 - nu, c, v); - /* Inu* */ + if (p == 2) + c = 1 + (a63 == 0) + ((a21 + a42 + a63) & 1); + else + c = 1 + (a63 == 0) + (((1 + a21 + a42 + a63) % 3) == 0) + + (((1 - a21 + a42 - a63) % 3) == 0); + return localred_result(nudelta - 4, -1, c, v); + /* I0* */ + case 2: /* compute nu */ + if (theroot) cumule(&v, &e, gun, stoi(theroot * p), gzero, gzero); + /* p | a1; p^2 | a2, a3; p^3 | a4; p^4 | a6 */ + nu = 1; + pk = stoi(p2); + p2k = stoi(p2 * p2); + for(;;) + { + be = aux2((GEN)e[3], p, pk); + ga = -aux2((GEN)e[5], p, p2k); + al = 1; + if (numroots2(al, be, ga, p, &theroot) == 2) break; + if (theroot) cumule(&v, &e, gun, gzero, gzero, mulsi(theroot,pk)); + pk1 = pk; + pk = mulsi(p, pk); + p2k = mulsi(p, p2k); nu++; + + al = a21; + be = aux2((GEN)e[4], p, pk); + ga = aux2((GEN)e[5], p, p2k); + if (numroots2(al, be, ga, p, &theroot) == 2) break; + if (theroot) cumule(&v, &e, gun, mulsi(theroot, pk1), gzero, gzero); + p2k = mulsi(p, p2k); nu++; + } + if (p == 2) + c = 4 - 2 * (ga & 1); + else + c = 3 + kross(be * be - al * ga, 3); + return localred_result(nudelta - 4 - nu, -4 - nu, c, v); + /* Inu* */ case 1: - if (theroot) cumule(&v, &e, gun, stoi(theroot * p), gzero, gzero); - /* p | a1; p^2 | a2, a3; p^3 | a4; p^4 | a6 */ - a32 = aux((GEN)e[3], p, 2); a64 = aux((GEN)e[5], p, 4); - if (numroots2(1, a32, -a64, p, &theroot) == 2) - { - if (p == 2) - c = 3 - 2 * a64; - else - c = 2 + kross(a32 * a32 + a64, 3); - return localreduction_result(av, nudelta - 6, -4, c, v); - } - /* IV* */ - if (theroot) cumule(&v, &e, gun, gzero, gzero, stoi(theroot*p*p)); - /* p | a1; p^2 | a2; p^3 | a3, a4; p^5 | a6 */ - p4 = sqri(p2); - if (!divise((GEN)e[4], p4)) - return localreduction_result(av, nudelta - 7, -3, 2, v); - /* III* */ - p6 = mulii(p4, p2); - if (!divise((GEN)e[5], p6)) - return localreduction_result(av, nudelta - 8, -2, 1, v); - /* II* */ - cumule(&v, &e, p1, gzero, gzero, gzero); /* non minimal, on repart - pour un tour */ - nudelta -= 12; + if (theroot) cumule(&v, &e, gun, stoi(theroot * p), gzero, gzero); + /* p | a1; p^2 | a2, a3; p^3 | a4; p^4 | a6 */ + a32 = aux((GEN)e[3], p, 2); + a64 = aux((GEN)e[5], p, 4); + if (numroots2(1, a32, -a64, p, &theroot) == 2) + { + if (p == 2) + c = 3 - 2 * a64; + else + c = 2 + kross(a32 * a32 + a64, 3); + return localred_result(nudelta - 6, -4, c, v); + } + /* IV* */ + if (theroot) cumule(&v, &e, gun, gzero, gzero, stoi(theroot*p*p)); + /* p | a1; p^2 | a2; p^3 | a3, a4; p^5 | a6 */ + p4 = p2 * p2; + if (smodis((GEN)e[4], p4)) + return localred_result(nudelta - 7, -3, 2, v); + /* III* */ + + if (smodis((GEN)e[5], p4 * p2)) + return localred_result(nudelta - 8, -2, 1, v); + /* II* */ + cumule(&v, &e, stoi(p), gzero, gzero, gzero); /* not minimal */ + nudelta -= 12; } } /* Not reached */ } +static GEN +localred(GEN e, GEN p, int minim) +{ + if (gcmpgs(p, 3) > 0) /* p != 2,3 */ + return localred_carac_p(e,p, minim); + else + { + GEN z = localred_carac_23(e,itos(p), minim); + return minim? (GEN)z[3]: z; + } +} + GEN -localreduction(GEN e, GEN p1) +localreduction(GEN e, GEN p) { + gpmem_t av = avma; checkell(e); if (typ(e[12]) != t_INT) err(talker,"not an integral curve in localreduction"); - if (gcmpgs(p1, 3) > 0) /* p different de 2 ou 3 */ - return localreduction_carac_not23(e,p1); - else - return localreduction_carac_23(e,p1); + return gerepileupto(av, localred(e, p, 0)); } -#if 0 -/* Calcul de toutes les fibres non elliptiques d'une courbe sur Z. - * Etant donne une courbe elliptique sous forme longue e, dont les coefficients - * sont entiers, renvoie une matrice dont les lignes sont de la forme - * [p, fp, kodp, cp]. Il y a une ligne par diviseur premier du discriminant. - */ -GEN -globaltatealgo(GEN e) +static GEN +ell_to_small(GEN E) { - long k, l,av; - GEN p1, p2, p3, p4, prims, result; + long i; + GEN e; + if (lg(E) <= 14) return E; + e = cgetg(14,t_VEC); + for (i = 1; i < 14; i++) e[i] = E[i]; + return e; +} +static GEN +ellintegralmodel(GEN e) +{ + GEN a = cgetg(6,t_VEC), v, prims, d, u; + long i, l, k; + checkell(e); - prims = decomp((GEN)e[12]); - l = lg(p1 = (GEN)prims[1]); - p2 = (GEN)prims[2]; - if ((long)prims == avma) cgiv(prims); - result = cgetg(5, t_MAT); - result[1] = (long)p1; - result[2] = (long)p2; - result[3] = (long)(p3 = cgetg(l, t_COL)); - for (k = 1; k < l; k++) p3[k] = lgeti(3); - result[4] = (long)(p4 = cgetg(l, t_COL)); - for (k = 1; k < l; k++) p4[k] = lgeti(3); - av = avma; + for (i = 1; i < 6; i++) + { + a[i] = e[i]; + switch(typ(a[i])) + { + case t_INT: case t_FRAC: case t_FRACN: break; + default: err(talker, "not a rational curve in ellintegralmodel"); + } + } + /* a = [a1, a2, a3, a4, a6] */ + d = denom(a); + prims = (GEN)auxdecomp(d, 0)[1]; /* partial factorization */ + if (is_pm1(d)) return NULL; + + l = lg(prims); u = gun; for (k = 1; k < l; k++) { - GEN q = localreduction(e, (GEN)p1[k]); - affii((GEN)q[1],(GEN)p2[k]); - affii((GEN)q[2],(GEN)p3[k]); - affii((GEN)q[4],(GEN)p4[k]); - avma = av; + GEN p = (GEN)prims[k]; + int n = 0, m; + for (i = 1; i < 6; i++) + if (!gcmp0((GEN)a[i])) + { + int r = (i == 5)? 6: i; /* a5 is missing */ + m = r * n + ggval((GEN)a[i], p); + while (m < 0) { n++; m += r; } + } + u = mulii(u, gpowgs(p, n)); } - return result; + v = init_ch(); v[1] = linv(u); return v; } -#endif -/* Algorithme de reduction d'une courbe sur Q a sa forme standard. Etant - * donne une courbe elliptique sous forme longue e, dont les coefficients - * sont rationnels, renvoie son [N, [u, r, s, t], c], ou N est le conducteur - * arithmetique de e, [u, r, s, t] est le changement de variables qui reduit - * e a sa forme minimale globale dans laquelle a1 et a3 valent 0 ou 1, et a2 - * vaut -1, 0 ou 1 et tel que u est un rationnel positif. Enfin c est le - * produit des nombres de Tamagawa locaux cp. - */ +static void +standard_model(GEN e, GEN *pv) +{ + GEN r, s, t; + s = gdiventgs((GEN)e[1], -2); + r = gdiventgs(gaddgs(gsub(gsub((GEN)e[2], gmul(s,(GEN)e[1])), gsqr(s)), 1), -3); + t = gdiventgs(ellLHS0(e,r), -2); + cumulev(pv, gun, r, s, t); +} + GEN -globalreduction(GEN e1) +ellminimalmodel(GEN E, GEN *ptv) { - long i, k, l, m, tetpil, av = avma; - GEN p1, c = gun, prims, result, N = gun, u = gun, r, s, t; - GEN v = cgetg(5, t_VEC), a = cgetg(7, t_VEC), e = cgetg(20, t_VEC); + gpmem_t av = avma; + GEN c4, c6, e, v, prims; + long l, k; - checkell(e1); - for (i = 1; i < 5; i++) a[i] = e1[i]; a[5] = zero; a[6] = e1[5]; - prims = decomp(denom(a)); - p1 = (GEN)prims[1]; l = lg(p1); + v = ellintegralmodel(E); + e = ell_to_small(E); + if (v) e = coordch(e, v); else v = init_ch(); + c4 = (GEN)e[10]; + c6 = (GEN)e[11]; + prims = (GEN)decomp(mppgcd(c4,c6))[1]; + l = lg(prims); for (k = 1; k < l; k++) { - int n = 0; - for (i = 1; i < 7; i++) - if (!gcmp0((GEN)a[i])) - { - m = i * n + ggval((GEN)a[i], (GEN)p1[k]); - while (m < 0) { n++; m += i; } - } - u = gmul(u, gpuigs((GEN)p1[k], n)); + GEN w = localred(e, (GEN)prims[k], 1); + if (!gcmp1((GEN)w[1])) + cumule(&v, &e, (GEN)w[1], (GEN)w[2], (GEN)w[3], (GEN)w[4]); } - v[1] = linv(u); v[2] = v[3] = v[4] = zero; - for (i = 1; i < 14; i++) e[i] = e1[i]; - for (; i < 20; i++) e[i] = zero; - if (!gcmp1(u)) e = coordch(e, v); - prims = decomp((GEN)e[12]); - l = lg(p1 = (GEN)prims[1]); - for (k = (signe(e[12]) < 0) + 1; k < l; k++) + standard_model(e, &v); + e = coordch(E, v); + if (ptv) { gerepileall(av, 2, &e, &v); *ptv = v; } + else e = gerepileupto(av, e); + return e; +} + +/* Reduction of a rational curve E to its standard minimal model + * (a1,a3 = 0 or 1, a2 = -1, 0 or 1). + * + * Return [N, [u,r,s,t], c], where + * N = arithmetic conductor of E + * c = product of the local Tamagawa numbers cp + * [u, r, s, t] = the change of variable reducing E to its minimal model, + * with u > 0 */ +GEN +globalreduction(GEN E) +{ + long k, l; + gpmem_t av = avma; + GEN c, prims, result, N, v, e; + + v = ellintegralmodel(E); + e = ell_to_small(E); + if (v) e = coordch(e, v); else v = init_ch(); + + prims = (GEN)decomp(absi((GEN)e[12]))[1]; + l = lg(prims); c = N = gun; + for (k = 1; k < l; k++) { - GEN q = localreduction(e, (GEN)p1[k]); - GEN v1 = (GEN)q[3]; - N = mulii(N, gpui((GEN)p1[k],(GEN)q[1],0)); + GEN p = (GEN)prims[k]; + GEN q = localreduction(e, p); + GEN w = (GEN)q[3]; + N = mulii(N, powgi(p, (GEN)q[1])); c = mulii(c, (GEN)q[4]); - if (!gcmp1((GEN)v1[1])) cumule1(&v, &e, v1); + if (!gcmp1((GEN)w[1])) + cumule(&v, &e, (GEN)w[1], (GEN)w[2], (GEN)w[3], (GEN)w[4]); } - s = gdiventgs((GEN)e[1], -2); - r = gdiventgs(gaddgs(gsub(gsub((GEN)e[2], gmul(s,(GEN)e[1])), gsqr(s)), 1), -3); - t = gdiventgs(ellLHS0(e,r), -2); - cumule(&v, &e, gun, r, s, t); - tetpil = avma; - result = cgetg(4, t_VEC); result[1] = lcopy(N); result[2] = lcopy(v); - result[3] = lcopy(c); - return gerepile(av, tetpil, result); + standard_model(e, &v); + result = cgetg(4, t_VEC); + result[1] = (long)N; + result[2] = (long)v; + result[3] = (long)c; return gerepilecopy(av, result); } -/* cumule les effets de plusieurs chgts de variable. On traite a part les cas - * particuliers frequents, tels que soit u = 1, soit r' = s' = t' = 0 - */ +/* accumulate the effects of variable changes [u,r,s,t] and [U,R,S,T]. + * Frequent special cases: (U = 1) or (r = s = t = 0) */ static void cumulev(GEN *vtotal, GEN u, GEN r, GEN s, GEN t) { - long av = avma, tetpil; - GEN temp, v = *vtotal, v3 = cgetg(5, t_VEC); - if (gcmp1((GEN)v[1])) + GEN U,R,S,T,UU,RR,SS,TT, v = *vtotal, w = cgetg(5, t_VEC); + gpmem_t av = avma; + + U = (GEN)v[1]; + R = (GEN)v[2]; + S = (GEN)v[3]; + T = (GEN)v[4]; + if (gcmp1(U)) { - v3[1] = lcopy(u); - v3[2] = ladd((GEN)v[2], r); - v3[3] = ladd((GEN)v[3], s); - av = avma; - temp = gadd((GEN)v[4], gmul((GEN)v[3], r)); - tetpil = avma; - v3[4] = lpile(av, tetpil, gadd(temp, t)); + UU = gcopy(u); + RR = gadd(R, r); + SS = gadd(S, s); av = avma; + TT = gerepileupto(av, gadd(T, gadd(t, gmul(S, r)))); } else if (gcmp0(r) && gcmp0(s) && gcmp0(t)) { - v3[1] = lmul((GEN)v[1], u); - v3[2] = lcopy((GEN)v[2]); - v3[3] = lcopy((GEN)v[3]); - v3[4] = lcopy((GEN)v[4]); + UU = gmul(U, u); + RR = gcopy(R); + SS = gcopy(S); + TT = gcopy(T); } - else /* cas general */ + else /* general case */ { - v3[1] = lmul((GEN)v[1], u); - temp = gsqr((GEN)v[1]); - v3[2] = ladd(gmul(temp, r), (GEN)v[2]); - v3[3] = ladd(gmul((GEN)v[1], s), (GEN)v[3]); - v3[4] = ladd((GEN)v[4], gmul(temp, gadd(gmul((GEN)v[1], t), gmul((GEN)v[3], r)))); - - v3 = gerepilecopy(av, v3); + GEN U2 = gsqr(U); + UU = gmul(U, u); + RR = gadd(R, gmul(U2, r)); + SS = gadd(S, gmul(U, s)); + TT = gadd(T, gmul(U2, gadd(gmul(U, t), gmul(S, r)))); + gerepileall(av, 4, &UU,&RR,&SS,&TT); } - *vtotal = v3; + w[1] = (long)UU; + w[2] = (long)RR; + w[3] = (long)SS; + w[4] = (long)TT; *vtotal = w; } static void cumule(GEN *vtotal, GEN *e, GEN u, GEN r, GEN s, GEN t) { - long av = avma, tetpil; - GEN v2 = cgetg(5, t_VEC); - v2[1] = (long)u; v2[2] = (long)r; v2[3] = (long)s; v2[4] = (long)t; - tetpil = avma; - *e = gerepile(av, tetpil, coordch(*e, v2)); + *e = _coordch(*e, u, r, s, t); cumulev(vtotal, u, r, s, t); } -static void -cumule1(GEN *vtotal, GEN *e, GEN v2) -{ - *e = coordch(*e, v2); - cumulev(vtotal, (GEN)v2[1], (GEN)v2[2], (GEN)v2[3], (GEN)v2[4]); -} - /********************************************************************/ /** **/ /** Parametrisation modulaire **/ @@ -2628,14 +2848,15 @@ GEN taniyama(GEN e) { GEN v,w,c,d,s1,s2,s3; - long n,m,av=avma,tetpil; + long n, m; + gpmem_t av=avma, tetpil; checkell(e); v = cgetg(precdl+3,t_SER); v[1] = evalsigne(1) | evalvalp(-2) | evalvarn(0); v[2] = un; c=gtoser(anell(e,precdl+1),0); setvalp(c,1); d=ginv(c); c=gsqr(d); - for (n=-3; n<=precdl-4; n++) + for (n=-3; n<=(long)precdl-4; n++) { if (n!=2) { @@ -2681,7 +2902,8 @@ GEN orderell(GEN e, GEN p) { GEN p1; - long av=avma,k; + long k; + gpmem_t av=avma; checkell(e); checkpt(p); k=typ(e[13]); @@ -2696,27 +2918,6 @@ orderell(GEN e, GEN p) avma=av; return gzero; } -/* one can do much better by factoring denom(D) (see ellglobalred) */ -static GEN -ellintegralmodel(GEN e) -{ - GEN a = cgetg(6,t_VEC), v; - long i; - - for (i=1; i<6; i++) - { - a[i]=e[i]; - switch(typ(a[i])) - { - case t_INT: case t_FRAC: case t_FRACN: break; - default: err(talker, "not a rational curve in ellintegralmodel"); - } - } - a = denom(a); if (gcmp1(a)) return NULL; - v = cgetg(5,t_VEC); - v[1]=linv(a); v[2]=v[3]=v[4]=zero; return v; -} - /* Using Lutz-Nagell */ /* p is a polynomial of degree exactly 3 with integral coefficients @@ -2771,9 +2972,9 @@ GEN torsellnagelllutz(GEN e) { GEN d,ld,pol,p1,lr,r,v,w,w2,w3; - long i,j,nlr,t,t2,k,k2,av=avma; + long i, j, nlr, t, t2, k, k2; + gpmem_t av=avma; - checkell(e); v = ellintegralmodel(e); if (v) e = coordch(e,v); pol = RHSpol(e); @@ -2854,20 +3055,21 @@ torsellnagelllutz(GEN e) /* Using Doud's algorithm */ -/* Input e and n, finds a bound for #Tor */ +/* finds a bound for #E_tor */ static long -torsbound(GEN e, long n) +torsbound(GEN e) { - long av = avma, m, b, c, d, prime = 2; + long m, b, c, prime = 2; + gpmem_t av = avma; byteptr p = diffptr; GEN D = (GEN)e[12]; - + long n = bit_accuracy(lgefint(D)) >> 3; + /* n = number of primes to try ~ 1 prime every 8 bits in D */ b = c = m = 0; p++; while (m -5) return NULL; @@ -2924,6 +3126,8 @@ best_in_cycle(GEN e, GEN p, long k) return (gsigne(d_ellLHS(e,p)) < 0)? invell(e,p): p; } +/* = E_tors, possibly NULL (= oo), p,q independant unless NULL + * order p = k, order q = 2 unless NULL */ static GEN tors(GEN e, long k, GEN p, GEN q, GEN v) { @@ -2939,7 +3143,6 @@ tors(GEN e, long k, GEN p, GEN q, GEN v) p = best_in_cycle(e,p,k); if (v) { - v[1] = linv((GEN)v[1]); p = pointch(p,v); q = pointch(q,v); } @@ -2954,10 +3157,7 @@ tors(GEN e, long k, GEN p, GEN q, GEN v) { p = best_in_cycle(e,p,k); if (v) - { - v[1] = linv((GEN)v[1]); p = pointch(p,v); - } r = cgetg(4,t_VEC); r[1] = lstoi(k); p1 = cgetg(2,t_VEC); p1[1] = r[1]; r[2] = (long)p1; p1 = cgetg(2,t_VEC); p1[1] = lcopy(p); @@ -2977,34 +3177,35 @@ tors(GEN e, long k, GEN p, GEN q, GEN v) GEN torselldoud(GEN e) { - long b,i,ord,av=avma,prec, k = 1; - GEN v,w1,w22,w1j,w12,p,tor1,tor2; + long b, i, ord, prec, k = 1; + gpmem_t av=avma; + GEN v,w,w1,w22,w1j,w12,p,tor1,tor2; checkbell(e); v = ellintegralmodel(e); if (v) e = coordch(e,v); - b = lgefint((GEN)e[12]) >> 1; /* b = size of sqrt(D) */ - prec = precision((GEN)e[15]); + b = DEFAULTPREC + ((lgefint((GEN)e[12])-2) >> 1); /* b >= size of sqrt(D) */ + w1 = (GEN)e[15]; + prec = precision(w1); if (prec < b) err(precer, "torselldoud"); - b = max(b, DEFAULTPREC); - if (b < prec) { prec = b; e = gprec_w(e, b); } - b = torsbound(e,3); + if (b < prec) { prec = b; e = gprec_w(e, b); w1 = (GEN)e[15]; } + b = torsbound(e); if (b==1) { avma=av; return tors(e,1,NULL,NULL, v); } + if (v) v[1] = linv((GEN)v[1]); w22 = gmul2n((GEN)e[16],-1); - w1 = (GEN)e[15]; if (b % 4) - { + { /* cyclic of order 1, p, 2p, p <= 5 */ p = NULL; for (i=10; i>1; i--) { if (b%i != 0) continue; w1j = gdivgs(w1,i); - p = torspnt(e,pointell(e,w1j,prec),i); + p = torspnt(e,w1j,i,prec); if (!p && i%2==0) { - p = torspnt(e,pointell(e,gadd(w22,w1j),prec),i); - if (!p) p = torspnt(e,pointell(e,gadd(w22,gmul2n(w1j,1)),prec),i); + p = torspnt(e,gadd(w22,w1j),i,prec); + if (!p) p = torspnt(e,gadd(w22,gmul2n(w1j,1)),i,prec); } if (p) { k = i; break; } } @@ -3012,60 +3213,67 @@ torselldoud(GEN e) } ord = 0; tor1 = tor2 = NULL; - w12 = gmul2n((GEN)e[15],-1); - if ((p = torspnt(e,pointell(e,w12,prec),2))) + w12 = gmul2n(w1,-1); + if ((p = torspnt(e,w12,2,prec))) { tor1 = p; ord++; } - if ((p = torspnt(e,pointell(e,w22,prec),2)) - || (!ord && (p = torspnt(e,pointell(e,gadd(w12,w22),prec),2)))) + w = w22; + if ((p = torspnt(e,w,2,prec))) { tor2 = p; ord += 2; } - + if (!ord) + { + w = gadd(w12,w22); + if ((p = torspnt(e,w,2,prec))) + { + tor2 = p; ord += 2; + } + } + p = NULL; switch(ord) { - case 0: + case 0: /* no point of order 2 */ for (i=9; i>1; i-=2) { if (b%i!=0) continue; - w1j=gdivgs((GEN)e[15],i); - p = torspnt(e,pointell(e,w1j,prec),i); + w1j = gdivgs(w1,i); + p = torspnt(e,w1j,i,prec); if (p) { k = i; break; } } break; - case 1: - p = NULL; + case 1: /* 1 point of order 2: w1 / 2 */ for (i=12; i>2; i-=2) { if (b%i!=0) continue; - w1j=gdivgs((GEN)e[15],i); - p = torspnt(e,pointell(e,w1j,prec),i); + w1j = gdivgs(w1,i); + p = torspnt(e,w1j,i,prec); if (!p && i%4==0) - p = torspnt(e,pointell(e,gadd(w22,w1j),prec),i); + p = torspnt(e,gadd(w22,w1j),i,prec); if (p) { k = i; break; } } if (!p) { p = tor1; k = 2; } break; - case 2: + case 2: /* 1 point of order 2: w = w2/2 or (w1+w2)/2 */ for (i=5; i>1; i-=2) { if (b%i!=0) continue; - w1j = gdivgs((GEN)e[15],i); - p = torspnt(e,pointell(e,gadd(w22,w1j),prec),i+i); + w1j = gdivgs(w1,i); + p = torspnt(e,gadd(w,w1j),2*i,prec); if (p) { k = 2*i; break; } } if (!p) { p = tor2; k = 2; } tor2 = NULL; break; - case 3: + case 3: /* 2 points of order 2: w1/2 and w2/2 */ for (i=8; i>2; i-=2) { if (b%(2*i)!=0) continue; - w1j=gdivgs((GEN)e[15],i); - p = torspnt(e,pointell(e,w1j,prec),i); + w1j = gdivgs(w1,i); + p = torspnt(e,w1j,i,prec); if (p) { k = i; break; } } if (!p) { p = tor1; k = 2; } @@ -3089,13 +3297,14 @@ elltors0(GEN e, long flag) /* par compatibilite */ GEN torsell(GEN e) {return torselldoud(e);} -/* LOCAL ROOT NUMBERS, D'APRES HALBERSTADT halberst@math.jussieu.fr */ +/* LOCAL ROOT NUMBERS, after HALBERSTADT */ -/* ici p=2 ou 3 */ +/* p = 2 or 3 */ static long neron(GEN e, GEN p, long* ptkod) { - long av=avma,kod,v4,v6,vd; + long kod, v4, v6, vd; + gpmem_t av=avma; GEN c4, c6, d, nv; nv=localreduction(e,p); @@ -3196,7 +3405,8 @@ neron(GEN e, GEN p, long* ptkod) static long ellrootno_2(GEN e) { - long n2,kod,u,v,x1,y1,d1,av=avma,v4,v6,w2; + long n2, kod, u, v, x1, y1, d1, v4, v6, w2; + gpmem_t av=avma; GEN p=gdeux,c4,c6,tmp,p6; n2=neron(e,p,&kod); c4=(GEN)e[10]; c6=(GEN)e[11]; p6=stoi(64); @@ -3283,7 +3493,8 @@ ellrootno_2(GEN e) static long ellrootno_3(GEN e) { - long n2,kod,u,v,d1,av=avma,r6,K4,K6,v4; + long n2, kod, u, v, d1, r6, K4, K6, v4; + gpmem_t av=avma; GEN p=stoi(3),c4,c6,tmp,p4; n2=neron(e,p,&kod); c4=(GEN)e[10]; c6=(GEN)e[11]; p4=stoi(81); @@ -3328,8 +3539,9 @@ ellrootno_3(GEN e) } } +/* assume p > 3 */ static long -ellrootno_not23(GEN e, GEN p, GEN ex) +ellrootno_p(GEN e, GEN p, GEN ex) { GEN j; long ep,z; @@ -3346,7 +3558,7 @@ ellrootno_not23(GEN e, GEN p, GEN ex) static long ellrootno_intern(GEN e, GEN p, GEN ex) { - if (cmpis(p,3) > 0) return ellrootno_not23(e,p,ex); + if (cmpis(p,3) > 0) return ellrootno_p(e,p,ex); switch(itos(p)) { case 3: return ellrootno_3(e); @@ -3378,7 +3590,7 @@ ellrootno_all(GEN e, GEN p, GEN* ptcond) exs=ggval(cond,p); if (!exs) return 1; } - if (cmpis(p,3)>0) return ellrootno_not23(e,p,stoi(exs)); + if (cmpis(p,3)>0) return ellrootno_p(e,p,stoi(exs)); switch(itos(p)) { case 3: return ellrootno_3(e); @@ -3393,7 +3605,8 @@ ellrootno_all(GEN e, GEN p, GEN* ptcond) long ellrootno(GEN e, GEN p) { - long av=avma,s; + long s; + gpmem_t av=avma; if (!p) p = gun; s=ellrootno_all(e, p, NULL); avma=av; return s;