home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Professional / OS2PRO194.ISO / os2 / progs / pari / pari_137 / src / elliptic.c < prev    next >
C/C++ Source or Header  |  1992-05-20  |  55KB  |  1,812 lines
  1. /********************************************************************/
  2. /********************************************************************/
  3. /**                                                                **/
  4. /**                    COURBES ELLIPTIQUES                         **/
  5. /**                                                                **/
  6. /**                     Copyright Babe Cool                        **/
  7. /**                                                                **/
  8. /********************************************************************/
  9. /********************************************************************/
  10.  
  11. #include "genpari.h"
  12.  
  13. static int aux(),aux2(),numroots2(),numroots3();
  14.  
  15. GEN smallinitell(x)
  16.      GEN x;
  17.  
  18. {
  19.   GEN y,b2,b4,b6,b8,d,j,a11,a13,a33,a64,b81,b22,c4,c6;
  20.   long i,av,tetpil;
  21.  
  22.   av=avma;y=cgetg(14,17);
  23.   if ((typ(x) != 17) || (lg(x) < 6)) err(elliper1);
  24.   for(i=1;i<=5;i++) y[i]=x[i];
  25.   b2=gadd(a11=gmul(y[1],y[1]),gmul2n(y[2],2));y[6]=(long)b2;
  26.   b4=gadd(a13=gmul(y[1],y[3]),gmul2n(y[4],1));y[7]=(long)b4;
  27.   b6=gadd(a33=gmul(y[3],y[3]),a64=gmul2n(y[5],2));y[8]=(long)b6;
  28.   b81=gadd(gadd(gmul(a11,y[5]),gmul(a64,y[2])),gmul(y[2],a33));
  29.   b8=gsub(b81,gmul(y[4],gadd(y[4],a13)));y[9]=(long)b8;
  30.   c4=gadd(b22=gmul(b2,b2),gmulsg(-24,b4));y[10]=(long)c4;
  31.   c6=gadd(gmul(b2,gsub(gmulsg(36,b4),b22)),gmulsg(-216,b6));y[11]=(long)c6;
  32.   b81=gadd(gmul(b22,b8),gmulsg(27,gmul(b6,b6)));
  33.   d=gsub(gmul(b4,gadd(gmulsg(9,gmul(b2,b6)),gmulsg(-8,gmul(b4,b4)))),b81);
  34.   y[12]=(long)d;j=gcmp0(d)?gzero:gdiv(gmul(gmul(c4,c4),c4),d);y[13]=(long)j;
  35.   tetpil=avma;return gerepile(av,tetpil,gcopy(y));
  36. }
  37.  
  38. static GEN initellcom(x,prec,hard)
  39.      GEN x;
  40.      long prec;
  41.      long hard;
  42.      
  43. {
  44.   GEN y,b2,b4,b6,b8,d,j,a11,a13,a33,a64,b81,b22,c4,c6,p1,p2,p,u,w,w2;
  45.   GEN aa0,bb0,aa1,bb1,r1,x0,x1,u2,q,pv,bmod0,bmod1,tau,pi2,pitemp;
  46.   long ty,i,av,tetpil,e0,e1,e4,ind,alpha,sw;
  47.  
  48.   av=avma;y=cgetg(20,17);
  49.   if ((typ(x) != 17) || (lg(x) < 6)) err(elliper1);
  50.   e0=BIGINT;
  51.   for(i=1;i<=5;i++) 
  52.     {
  53.       q=(GEN)x[i];y[i]=(long)q;
  54.       if(typ(q)==7) 
  55.     {
  56.       e0=min(e0,signe(q[4])?precp(q)+valp(q):valp(q));
  57.       p=(GEN)q[2];
  58.     }
  59.     }
  60.   if(e0<BIGINT)
  61.     {
  62.       q=ggrando(p,e0);for(i=1;i<=5;i++) y[i]=ladd(q,x[i]);
  63.     }
  64.   b2=gadd(a11=gmul(y[1],y[1]),gmul2n(y[2],2));y[6]=(long)b2;
  65.   b4=gadd(a13=gmul(y[1],y[3]),gmul2n(y[4],1));y[7]=(long)b4;
  66.   b6=gadd(a33=gmul(y[3],y[3]),a64=gmul2n(y[5],2));y[8]=(long)b6;
  67.   b81=gadd(gadd(gmul(a11,y[5]),gmul(a64,y[2])),gmul(y[2],a33));
  68.   b8=gsub(b81,gmul(y[4],gadd(y[4],a13)));y[9]=(long)b8;
  69.   c4=gadd(b22=gmul(b2,b2),gmulsg(-24,b4));y[10]=(long)c4;
  70.   c6=gadd(gmul(b2,gsub(gmulsg(36,b4),b22)),gmulsg(-216,b6));y[11]=(long)c6;
  71.   b81=gadd(gmul(b22,b8),gmulsg(27,gmul(b6,b6)));
  72.   d=gsub(gmul(b4,gadd(gmulsg(9,gmul(b2,b6)),gmulsg(-8,gmul(b4,b4)))),b81);
  73.   y[12]=(long)d;j=gdiv(gmul(gmul(c4,c4),c4),d);y[13]=(long)j;
  74.   if(gcmp0(d)) err(initeler2);
  75.   ty=typ(d);
  76.   if(prec&&(ty<=8)&&(ty!=3)&&(ty!=7))
  77.     {
  78.       p1=cgetg(6,10);p1[1]=0x1000006;p1[2]=(long)b6;
  79.       p1[3]=lmul2n(b4,1);p1[4]=(long)b2;p1[5]=lstoi(4);
  80.       if(hard) p1=rootslong(p1,prec);else p1=roots(p1,prec);
  81.       if(gsigne(d)>0) 
  82.         {
  83.           p1=greal(p1);
  84.           ind=1;e4=p1[1];for(i=2;i<=3;i++) 
  85.         if(gcmp(p1[i],e4)>0) {ind=i;e4=p1[i];}
  86.           p1[ind]=p1[1];p1[1]=e4;
  87.           if(gcmp(p1[2],p1[3])<0) {e4=p1[3];p1[3]=p1[2];p1[2]=e4;}
  88.         }
  89.       else p1[1]=lreal(p1[1]);
  90.       y[14]=(long)p1;e1=p1[1];p2=gadd(gmulsg(3,e1),gmul2n(b2,-2));
  91.       w=gsqrt(gmul2n(gadd(b4,gmul(e1,gadd(b2,gmulsg(6,e1)))),1),prec);
  92.       if(gsigne(p2)>0) w=gneg(w);
  93.       aa1=gmul2n(gsub(w,p2),-2);sw=signe(w);
  94.       bb1=gmul2n(w,-1);r1=gsub(aa1,bb1);x1=gmul2n(r1,-2);
  95.       do
  96.         {
  97.           aa0=aa1;bb0=bb1;x0=x1;
  98.           bb1=gsqrt(gmul(aa0,bb0),prec);setsigne(bb1,sw);
  99.           aa1=gmul2n(gadd(gadd(aa0,bb0),gmul2n(bb1,1)),-2);
  100.           r1=gsub(aa1,bb1);
  101.       x1=gmul(x0,gsqr(gmul2n(gaddsg(1,gsqrt(gdiv(gadd(x0,r1),x0),prec)),-1)));
  102.         }
  103.       while(gexpo(r1)>=gexpo(bb1)-((prec-2)<<5)+5);
  104.       u2=ginv(gmul2n(aa1,2));
  105.       w=gaddsg(1,ginv(gmul2n(gmul(u2,x1),1)));
  106.       if(gsigne(greal(w))>0) q=ginv(gadd(w,gsqrt(gaddgs(gmul(w,w),-1),prec)));
  107.       else q=gsub(w,gsqrt(gaddgs(gmul(w,w),-1),prec));
  108.       pitemp=mppi(prec);pi2=gmul2n(pitemp,1);if(gexpo(q)>=0) q=ginv(q);
  109.       tau=gmul(gdiv(glog(q,prec),pi2),gneg(gi));
  110.       y[19]=lmul(gmul(gsqr(pi2),gabs(u2,prec)),gimag(tau));
  111.       u=gmul(pi2,gsqrt(gneg(u2),prec));w2=gmul(tau,u);
  112.       if(sw<0) y[15]=(long)u;
  113.       else
  114.         {
  115.           y[15]=lmul2n(gabs(w2[1],prec),1);
  116.           q=gexp(gmul(gmul(pi2,gi),gdiv(w2,y[15])),prec);
  117.         }
  118.       y[16]=(long)w2;
  119.       p1=gdiv(gmul(pitemp,pitemp),gmulsg(6,y[15]));q=gsqrt(q,prec);
  120.       y[17]=lmul(p1,gdiv(thetanullk(q,3,prec),thetanullk(q,1,prec)));
  121.       y[18]=ldiv(gsub(gmul(y[17],y[16]),gmul(gi,pitemp)),y[15]);
  122.     }
  123.   else 
  124.     if(ty==7)
  125.       {
  126.         if(valp(j)>=0) err(initeler1);
  127.         p=(GEN)d[2];alpha=valp(c4)>>1;
  128.         setvalp(c4,0);setvalp(c6,0);e1=ldivgs(gdiv(c6,c4),6);
  129.         c4=gdivgs(c4,48);c6=gdivgs(c6,864);
  130.         do
  131.           {
  132.             e0=e1;p2=gmul(e0,e0);
  133.             e1=ldiv(gadd(gmul2n(gmul(e0,p2),1),c6),gsub(gmulsg(3,p2),c4));
  134.           }
  135.         while(!gegal(e0,e1));
  136.         setvalp(e1,valp(e1)+alpha);e1=lsub(e1,gdivgs(b2,12));
  137.         w=gsqrt(gmul2n(gadd(b4,gmul(e1,gadd(b2,gmulsg(6,e1)))),1));
  138.         p1=gaddgs(gdiv(gmulsg(3,e0),w),1);
  139.         if(valp(p1)<=0) w=gneg(w);
  140.         y[18]=(long)w;
  141.         aa1=gmul2n(gsub(w,gadd(gmulsg(3,e1),gmul2n(b2,-2))),-2);
  142.         bb1=gmul2n(w,-1);r1=gsub(aa1,bb1);x1=gmul2n(r1,-2);
  143.     if(gegal(p,deux)) {pv=stoi(4);err(impl,"initell for 2-adic numbers");}
  144.     else pv=p;
  145.         bmod0=modii(bb1[4],pv);
  146.         do
  147.           {
  148.             aa0=aa1;bb0=bb1;x0=x1;
  149.             bb1=gsqrt(gmul(aa0,bb0));bmod1=modii(bb1[4],pv);
  150.             if(!gegal(bmod1,bmod0)) bb1=gneg(bb1);
  151.             aa1=gmul2n(gadd(gadd(aa0,bb0),gmul2n(bb1,1)),-2);
  152.             r1=gsub(aa1,bb1);
  153.             p1=gsqrt(gdiv(gadd(x0,r1),x0));
  154.             if(!gegal(modii(p1[4],pv),un)) p1=gneg(p1);
  155.             x1=gmul(x0,gsqr(gmul2n(gaddsg(1,p1),-1)));
  156.           }
  157.         while(!gcmp0(r1));
  158.         u2=ginv(gmul2n(aa1,2));
  159.         w=gaddsg(1,ginv(gmul2n(gmul(u2,x1),1)));
  160.         q=ginv(gadd(w,gsqrt(gaddgs(gmul(w,w),-1))));
  161.         if(valp(q)<0) q=ginv(q);
  162.         y[15]=(long)u2;
  163.         y[16]=((kronecker(u2[4],p)>0)&&(!(valp(u2)&1)))?lsqrt(u2):zero;
  164.         y[17]=(long)q;p1=cgetg(2,17);
  165.         y[14]=(long)p1;p1[1]=e1;y[19]=zero;
  166.       }
  167.     else  {y[14]=y[15]=y[16]=y[17]=y[18]=y[19]=zero;}
  168.   tetpil=avma;return gerepile(av,tetpil,gcopy(y));
  169. }
  170.  
  171. GEN initell2(x,prec)
  172.      GEN x;
  173.      long prec;
  174.  
  175. {
  176.   return initellcom(x, prec, 0);
  177. }
  178.  
  179. GEN initell(x,prec)
  180.      GEN x;
  181.      long prec;
  182.  
  183. {
  184.   return initellcom(x, prec, 1);
  185. }
  186.  
  187. GEN coordch(e,ch)
  188.      GEN e,ch;
  189. {
  190.   GEN y,p1,p2,v,v2,v3,v4,v6,r,s,t,u;
  191.   long av,tetpil,i,lx=lg(e);
  192.   
  193.   if ((typ(e) != 17) || (lx < 14)) err(elliper1);
  194.   u=(GEN)ch[1];r=(GEN)ch[2];s=(GEN)ch[3];t=(GEN)ch[4];
  195.   av=avma;y=cgetg(lx,17);v=ginv(u);
  196.   y[1]=lmul(v,gadd(e[1],gmul2n(s,1)));v2=gmul(v,v);
  197.   y[2]=lmul(v2,gsub(gadd(e[2],gmulsg(3,r)),gmul(s,gadd(e[1],s))));
  198.   v3=gmul(v,v2);
  199.   y[3]=lmul(v3,p1=gadd(gadd(gmul(r,e[1]),gmul2n(t,1)),e[3]));
  200.   v4=gmul(v,v3);
  201.   p1=gsub(e[4],gadd(gmul(t,e[1]),gmul(s,p1)));
  202.   y[4]=lmul(v4,gadd(p1,gmul(r,gadd(gmul2n(e[2],1),gmulsg(3,r)))));
  203.   p1=gadd(e[5],gmul(r,gadd(e[4],gmul(r,gadd(e[2],r)))));
  204.   p2=gmul(t,gadd(gadd(gmul(r,e[1]),t),e[3]));v6=gmul(v2,v4);
  205.   y[5]=lmul(v6,gsub(p1,p2));
  206.   y[6]=lmul(v2,gadd(e[6],gmulsg(12,r)));
  207.   y[7]=lmul(v4,gadd(e[7],gmul(r,gadd(e[6],gmulsg(6,r)))));
  208.  
  209. y[8]=lmul(v6,gadd(e[8],gmul(r,gadd(gmul2n(e[7],1),gmul(r,gadd(e[6],gmul2n(r,2)
  210. ))))));
  211.   p1=gadd(gmulsg(3,e[7]),gmul(r,gadd(e[6],gmulsg(3,r))));
  212.   y[9]=lmul(gmul(v6,v2),gadd(e[9],gmul(r,gadd(gmulsg(3,e[8]),gmul(r,p1)))));
  213.   y[10]=lmul(v4,e[10]);y[11]=lmul(v6,e[11]);
  214.   y[12]=lmul(gmul(v6,v6),e[12]);y[13]=lcopy(e[13]);
  215.   if(lx==14) {tetpil=avma;return gerepile(av,tetpil,gcopy(y));}
  216.   p1=(GEN)e[14];
  217.   if (gcmp0(p1))
  218.   {
  219.     y[14] = y[15] = y[16] = y[17] = y[18] = y[19] = zero;
  220.   }
  221.   else
  222.   {
  223.     p2=cgetg(4,18);
  224.     for(i=1;i<=3;i++) p2[i]=lmul(v2,gsub(p1[i],r));
  225.     y[14]=(long)p2;
  226.     y[15]=lmul(u,e[15]);if(typ(e[1])==7) y[15]=lmul(u,y[15]);
  227.     y[16]=lmul(u,e[16]);
  228.     y[17]=ldiv(e[17],u);y[18]=ldiv(e[18],u);
  229.     y[19]=lmul(gmul(u,u),e[19]);
  230.   }
  231.   tetpil=avma;return gerepile(av,tetpil,gcopy(y));
  232. }
  233.  
  234. GEN pointch(x,ch)
  235.      GEN x,ch;
  236. {
  237.   GEN y,p1,p2,p3,v,v2,v3,mor,r,s,t,u;
  238.   long av,tetpil,tx,lx=lg(x),i;
  239.  
  240.   if ((typ(x) != 17) || (typ(ch) != 17)) err(elliper1);
  241.   if (lx < 2) return gcopy(x);
  242.   av=avma;u=(GEN)ch[1];r=(GEN)ch[2];s=(GEN)ch[3];t=(GEN)ch[4];
  243.   tx=typ(x[1]);v=ginv(u);v2=gmul(v,v);v3=gmul(v,v2);
  244.   mor=gneg(r);
  245.   if(tx>=17)
  246.     {
  247.       y=cgetg(lx,tx);
  248.       for(i=1;i<lx;i++)
  249.         {
  250.           p2=(GEN)x[i];
  251.           if(lg(p2)<3) y[i]=(long)p2;
  252.           else
  253.             {
  254.               p1=cgetg(3,17);
  255.               p1[1]=lmul(v2,p3=gadd(p2[1],mor));
  256.               p1[2]=lmul(v3,gsub(p2[2],gadd(gmul(s,p3),t)));
  257.               y[i]=(long)p1;
  258.             }
  259.         }
  260.     }
  261.   else
  262.     if(lg(x)<3) y=x;
  263.     else
  264.       {
  265.     y=cgetg(3,17);
  266.     y[1]=lmul(v2,p3=gadd(x[1],mor));
  267.     y[2]=lmul(v3,gsub(x[2],gadd(gmul(s,p3),t)));
  268.       }
  269.   tetpil=avma;return gerepile(av,tetpil,gcopy(y));
  270. }
  271.  
  272. int oncurve(e,z)
  273.      GEN e,z;
  274. {
  275.   GEN p1,p2,x,y;
  276.   long av=avma,f;
  277.  
  278.   if ((typ(e) != 17) || (lg(e) < 6) || (typ(z) != 17)) err(elliper1);
  279.   if(lg(z)<3) return 1; /* cas du point a l'infini */
  280.   x=(GEN)z[1];y=(GEN)z[2];
  281.   p1=gmul(y,gadd(y,gadd(e[3],gmul(e[1],x))));
  282.   p2=gadd(e[5],gmul(x,gadd(e[4],gmul(x,gadd(e[2],x)))));
  283.   p1=gsub(p1,p2);
  284.   f=precision(p1);
  285.   if(!f) {f=gcmp0(p1);avma=av;return f;}
  286.   f=gcmp0(p1)||(gexpo(p1)< -((f-2)<<5)+20);avma=av;return f;
  287. }
  288.  
  289. GEN hells(e,x,prec)
  290.      GEN e,x;
  291.      long prec;
  292. {
  293.   GEN w,z,t,mu,unreel;
  294.   long av=avma,tetpil,n,lim;
  295.  
  296.   if(!oncurve(e,x)) err(heller1);
  297.   affsr(1,unreel=cgetr(prec));
  298.   t=gdiv(unreel,x[1]);mu=gmul2n(glog(numer(x[1]),prec),-1);
  299.   n=0;lim=(prec-2)<<4;
  300.   do
  301.     {
  302.      
  303. w=gmul(t,gaddsg(4,gmul(t,gadd(e[6],gmul(t,gadd(gmul2n(e[7],1),gmul(t,e[8])))))
  304. ));
  305.      
  306. z=gsub(gun,gmul(gmul(t,t),gadd(e[7],gmul(t,gadd(gmul2n(e[8],1),gmul(t,e[9]))))
  307. ));
  308.       mu=gadd(mu,gmul2n(glog(z,prec),-((n<<1)+3)));t=gdiv(w,z);n++;
  309.     }
  310.   while(n<=(lim+5));
  311.   tetpil=avma;return gerepile(av,tetpil,gcopy(mu));
  312. }
  313.  
  314. GEN hell2(e,x,prec)
  315.      GEN e,x;
  316.      long prec;
  317. {
  318.   GEN ep,e3,ro,p1,p2,mu,d,xp;
  319.   long av=avma,tetpil,lx,lc,i,j,tx;
  320.  
  321.   if((typ(e)!=17)||(lg(e)<20)) err(elliper1);
  322.   if(!oncurve(e,x)) err(heller1);
  323.   d=(GEN)e[12];ro=(GEN)e[14];e3=(gsigne(d)<0)?(GEN)ro[1]:(GEN)ro[3];
  324.   p1=cgetg(5,17);p1[1]=un;p1[2]=laddgs(gfloor(e3),-1);
  325.   p1[3]=p1[4]=zero;ep=coordch(e,p1);xp=pointch(x,p1);
  326.   if(typ(x[1])<17)
  327.     {
  328.       if(lg(x)<3) return gzero;
  329.       tetpil=avma;return gerepile(av,tetpil,hells(ep,xp,prec));
  330.     }
  331.   else
  332.     {
  333.       lx=lg(x);tetpil=avma;mu=cgetg(lx,tx=typ(x));
  334.       if(tx<19) for(i=1;i<lx;i++) mu[i]=(long)hells(ep,xp[i],prec);
  335.       else
  336.         {
  337.           lc=lg(x[1]);
  338.           for(i=1;i<lx;i++)
  339.             {
  340.               p1=cgetg(lc,18);mu[i]=(long)p1;p2=(GEN)xp[i];
  341.               for(j=1;j<lc;j++) p1[j]=(long)hells(ep,p2[j],prec);
  342.             }
  343.         }
  344.       return gerepile(av,tetpil,mu);
  345.     }
  346. }
  347.  
  348. GEN addell(e,z1,z2)
  349.      GEN e,z1,z2;
  350.   
  351. {
  352.   GEN y,p1,p2,x1,x2,x3,y1,y2,y3,al;
  353.   long av=avma,tetpil;
  354.  
  355.   if((typ(e) != 17) || (typ(z1) != 17) || (typ(z2) != 17)) err(elliper1);
  356.   if(lg(z1)<3) return gcopy(z2);
  357.   if(lg(z2)<3) return gcopy(z1);
  358.   x1=(GEN)z1[1];x2=(GEN)z2[1];y1=(GEN)z1[2];y2=(GEN)z2[2];
  359.   if(gegal(x1,x2))
  360.     if(!gegal(y1,y2)) {y=cgetg(2,17);y[1]=zero;return y;}
  361.     else
  362.       {
  363.        
  364. p1=gadd(gsub(e[4],gmul(e[1],y1)),gmul(x1,gadd(gmul2n(e[2],1),gmulsg(3,x1))));
  365.         p2=gadd(e[3],gadd(gmul(e[1],x1),gmul2n(y1,1)));
  366.         if(gcmp0(p2)) {avma=av;y=cgetg(2,17);y[1]=zero;return y;}
  367.       }
  368.   else {p1=gsub(y2,y1);p2=gsub(x2,x1);}
  369.   al=gdiv(p1,p2);
  370.   x3=gsub(gmul(al,gadd(al,e[1])),gadd(gadd(x1,x2),e[2]));
  371.   y3=gneg(gadd(gadd(gadd(e[3],y1),gmul(x3,e[1])),gmul(al,gsub(x3,x1))));
  372.   tetpil=avma;y=cgetg(3,17);y[1]=lcopy(x3);y[2]=lcopy(y3);
  373.   return gerepile(av,tetpil,y);
  374. }
  375.  
  376. GEN subell(e,z1,z2)
  377.      GEN e,z1,z2;
  378.  
  379. {
  380.   GEN zp;
  381.   long av=avma,tetpil;
  382.  
  383.   if((typ(e) != 17) || (typ(z2) != 17)) err(elliper1);
  384.   if(lg(z2)<3) return gcopy(z1);
  385.   zp=cgetg(3,17);zp[1]=z2[1];
  386.   zp[2]=lneg(gadd(gadd(z2[2],e[3]),gmul(e[1],z2[1])));
  387.   tetpil=avma;return gerepile(av,tetpil,addell(e,z1,zp));
  388. }
  389.  
  390. GEN ordell(e,x,prec)
  391.      GEN e,x;
  392.      long prec;
  393. {
  394.   GEN p1,p2,p3,d,y,pn,pd;
  395.   long av=avma,tetpil,td,i,lx,tx=typ(x);
  396.   
  397.   if((typ(e) != 17) || (lg(e) < 6)) err(elliper1);
  398.   if(tx<17)
  399.     {
  400.       p1=gadd(e[5],gmul(x,gadd(e[4],gmul(x,gadd(e[2],x)))));
  401.       p2=gadd(e[3],gmul(e[1],x));d=gadd(gmul(p2,p2),gmul2n(p1,2));
  402.       if((td=typ(d))==1)
  403.         {
  404.           if(signe(d)<0) {avma=av;return cgetg(1,17);}
  405.           if(!carrecomplet(d,&p3)) {avma=av;return cgetg(1,17);} 
  406.           p1=gsub(p3,p2);tetpil=avma;
  407.           if(gsigne(d)) {y=cgetg(3,17);y[1]=lmul2n(p1,-1);y[2]=lsub(y[1],p3);}
  408.           else {y=cgetg(2,17);y[1]=lmul2n(p1,-1);}
  409.           return gerepile(av,tetpil,y);
  410.         }
  411.       if((td==4)||(td==5))
  412.     {
  413.       pn=(GEN)d[1];pd=(GEN)d[2];
  414.       if(!carrecomplet(mulii(pn,pd),&p3)) {avma=av;return cgetg(1,17);} 
  415.           p1=gsub(p3=gdiv(p3,pd),p2);tetpil=avma;
  416.           if(signe(pn)) {y=cgetg(3,17);y[1]=lmul2n(p1,-1);y[2]=lsub(y[1],p3);}
  417.           else {y=cgetg(2,17);y[1]=lmul2n(p1,-1);}
  418.           return gerepile(av,tetpil,y);
  419.         }
  420.       if((td==3)&&gegal(d[1],gdeux))
  421.     {
  422.       if(gcmp0(p2)) 
  423.         {
  424.           avma=av;y=cgetg(2,17);p3=cgetg(3,3);y[1]=(long)p3;
  425.           p3[1]=deux;p3[2]=gcmp0(p1)?zero:un;
  426.         }
  427.       else
  428.         if(gcmp0(p1))
  429.           {
  430.         avma=av;y=cgetg(3,17);p3=cgetg(3,3);y[1]=(long)p3;
  431.         p3[1]=deux;p3[2]=zero;p3=cgetg(3,3);y[2]=(long)p3;
  432.         p3[1]=deux;p3[2]=un;
  433.           }
  434.         else {avma=av;y=cgetg(1,17);}
  435.       return y;
  436.     }
  437.       else
  438.     {
  439.       if((td==3)&&(kronecker(d[2],d[1])== -1))
  440.         {avma=av;y=cgetg(1,17);return y;}
  441.       p3=gsqrt(d,prec);p1=gsub(p3,p2);tetpil=avma;
  442.       if(!gcmp0(d))
  443.         {
  444.           y=cgetg(3,17);y[1]=lmul2n(p1,-1);y[2]=lsub(y[1],p3);
  445.         }
  446.       else {y=cgetg(2,17);y[1]=lmul2n(p1,-1);}
  447.       return gerepile(av,tetpil,y);
  448.     }
  449.     }
  450.   else
  451.     {
  452.       lx=lg(x);y=cgetg(lx,tx);
  453.       for(i=1;i<lx;i++) y[i]=(long)ordell(e,x[i],prec);
  454.       return y;
  455.     }
  456. }
  457.  
  458. GEN powell(e,n,z)
  459.      GEN e,n,z;
  460.  
  461. {
  462.   GEN y,zp;
  463.   long s=signe(n),av=avma,i,j,tetpil;
  464.   unsigned long m;
  465.  
  466.   if(typ(e) != 17) err(elliper1);
  467.   if(!s) {y=cgetg(2,17);y[1]=zero;return y;}
  468.   if(lg(z)<3) return gcopy(z);
  469.   if(s<0) 
  470.     {
  471.       n=gneg(n);zp=cgetg(3,17);zp[1]=z[1];
  472.       zp[2]=lneg(gadd(gadd(z[2],e[3]),gmul(e[1],z[1])));
  473.     }
  474.   else zp=z;
  475.   if(gcmp1(n)) {tetpil=avma;return gerepile(av,tetpil,gcopy(zp));}
  476.   y=cgetg(2,17);y[1]=zero;
  477.   for (i=lgef(n)-1;i>2;i--)
  478.     for (m=n[i],j=0;j<32;j++,m>>=1)
  479.       {
  480.         if (m&1) y=addell(e,y,zp);
  481.         zp=addell(e,zp,zp);
  482.       }
  483.   for (m=n[2];m>1;m>>=1)
  484.     {
  485.       if (m&1) y=addell(e,y,zp);
  486.       zp=addell(e,zp,zp);
  487.     }
  488.   tetpil=avma;y=addell(e,y,zp);
  489.   return gerepile(av,tetpil,y);
  490. }
  491.  
  492. GEN matell(e,x,prec)
  493.      GEN e,x;
  494.      long prec;
  495. {
  496.   GEN y,p1,p2,pdiag;
  497.   long av=avma,tetpil,lx=lg(x),i,j;
  498.  
  499.   if(typ(x) != 17) err(elliper1);
  500.   lx=lg(x);y=cgetg(lx,19);pdiag=cgetg(lx,18);
  501.   for(i=1;i<lx;i++) {pdiag[i]=(long)ghell(e,x[i],prec);y[i]=lgetg(lx,18);}
  502.   for(i=1;i<lx;i++)
  503.     {
  504.       p1=(GEN)y[i];p1[i]=lmul2n(pdiag[i],1);
  505.       for(j=i+1;j<lx;j++)
  506.         {
  507.       p2=gsub(ghell(e,addell(e,x[i],x[j]),prec),gadd(pdiag[i],pdiag[j]));
  508.           p1[j]=(long)p2;((GEN)y[j])[i]=(long)p2;
  509.         }
  510.     }
  511.   tetpil=avma;return gerepile(av,tetpil,gcopy(y));
  512. }
  513.      
  514. GEN zell(e,z,prec)
  515.      GEN e,z;
  516.      long prec;
  517. {
  518.   long av=avma,tetpil,ty,sw,fl;
  519.   GEN t,u,w,p1,p2,b2,b4,r0,r1,aa1,aa0,bb1,bb0,x0,x1,c0,e1,delta;
  520.   GEN bmod0,bmod1,p,u2;
  521.  
  522.   if((typ(e)!=17)||(lg(e)<20)) err(elliper1);
  523.   if(!oncurve(e,z)) err(heller1);
  524.   ty=typ(e[12]);b2=(GEN)e[6];b4=(GEN)e[7];
  525.   if(lg(z)<3) return (ty==7)?gun:gzero;
  526.   e1=(GEN)((GEN)e[14])[1];
  527.   w=gsqrt(gmul2n(gadd(b4,gmul(e1,gadd(b2,gmulsg(6,e1)))),1),prec);
  528.   if((ty<=8)&&(ty!=3)&&(ty!=7))
  529.     {
  530.       p2=gadd(gmulsg(3,e1),r0=gmul2n(b2,-2));
  531.       if(gsigne(p2)>0) w=gneg(w);
  532.       aa1=gmul2n(gsub(w,p2),-2);sw=signe(w);r0=gadd(e1,r0);
  533.       bb1=gmul2n(w,-1);r1=gsub(aa1,bb1);
  534.       c0=gadd(z[1],gmul2n(r0,-1));
  535.       if(gsigne(c0))
  536.     {
  537.       delta=gdiv(gmul(aa1,r1),gsqr(c0));
  538.       x0=gmul2n(gmul(c0,gaddsg(1,gsqrt(gsubsg(1,gmul2n(delta,2)),prec))),-1);
  539.     }
  540.       else x0=gsqrt(gneg(gmul(aa1,r1)),prec);
  541.       x1=gmul(x0,gsqr(gmul2n(gaddsg(1,gsqrt(gdiv(gadd(x0,r1),x0),prec)),-1)));
  542.       for(fl=0;;)
  543.         {
  544.           aa0=aa1;bb0=bb1;x0=x1;r0=r1;
  545.           bb1=gsqrt(gmul(aa0,bb0),prec);setsigne(bb1,sw);
  546.           aa1=gmul2n(gadd(gadd(aa0,bb0),gmul2n(bb1,1)),-2);
  547.           r1=gsub(aa1,bb1);
  548.       x1=gmul(x0,gsqr(gmul2n(gaddsg(1,gsqrt(gdiv(gadd(x0,r1),x0),prec)),-1)));
  549.       if(gexpo(gsub(x1,x0))<gexpo(x1)-((prec-2)<<5)+5)
  550.         if (fl) break; else fl = 1;
  551.       else fl = 0;
  552.         }
  553.       t=gaddsg(1,u=gdiv(x1,aa1));
  554.       if(gcmp0(t)||(gexpo(t)<-((prec-2)<<5)+5)) t=gneg(gun);
  555.       else t=gdiv(u,gsqr(gaddsg(1,gsqrt(t,prec))));
  556.       u=gsqrt(ginv(gmul2n(aa1,2)),prec);t=glog(t,prec);t=gmul(t,u);
  557. /* les deux lignes suivantes ont ete rajoutees au pif. A verifier */
  558.       if(gsigne(c0)*gsigne(gadd(e[3],gadd(gmul(e[1],z[1]),gmul2n(z[2],1))))<0)
  559.     t=gneg(t);
  560.       p1=gsub(p2=gdiv(gimag(t),((GEN)e[16])[2]),gmul2n(gun,-2));
  561.       if(gcmp(gabs(p1,prec),ghalf)>=0)
  562.         {
  563.           t=gsub(t,gmul(e[16],gfloor(gadd(p2,lisexpr("0.1")))));
  564.         }
  565.       if(signe(greal(t))<0) t=gadd(t,e[15]);
  566.       tetpil=avma;return gerepile(av,tetpil,gcopy(t));
  567.     }
  568.   else 
  569.     if(ty==7)
  570.       {
  571.         w=(GEN)e[18];r0=gadd(e1,gmul2n(b2,-2));p=(GEN)((GEN)e[12])[2];
  572.         aa1=gmul2n(gsub(w,gadd(gmulsg(3,e1),gmul2n(b2,-2))),-2);
  573.         bb1=gmul2n(w,-1);r1=gsub(aa1,bb1);
  574.         c0=gadd(z[1],gmul2n(r0,-1));delta=gdiv(gmul(aa1,r1),gsqr(c0));
  575.     x0=gmul2n(gmul(c0,gaddsg(1,gsqrt(gsubsg(1,gmul2n(delta,2)),prec))),-1);
  576.         bmod0=modii(bb1[4],p);
  577.     x1=gmul(x0,gsqr(gmul2n(gaddsg(1,gsqrt(gdiv(gadd(x0,r1),x0),prec)),-1)));
  578.         do
  579.           {
  580.             aa0=aa1;bb0=bb1;x0=x1;r0=r1;
  581.             bb1=gsqrt(gmul(aa0,bb0));bmod1=modii(bb1[4],p);
  582.             if(!gegal(bmod1,bmod0)) bb1=gneg(bb1);
  583.             aa1=gmul2n(gadd(gadd(aa0,bb0),gmul2n(bb1,1)),-2);
  584.             r1=gsub(aa1,bb1);
  585.             p1=gsqrt(gdiv(gadd(x0,r1),x0));
  586.             if(!gegal(modii(p1[4],p),un)) p1=gneg(p1);
  587.             x1=gmul(x0,gsqr(gmul2n(gaddsg(1,p1),-1)));
  588.           }
  589.         while(!gcmp0(r1));
  590.         u2=ginv(gmul2n(aa1,2));
  591.         if(!gcmp0(e[16]))
  592.           {
  593.             t=gsqrt(gaddsg(1,gdiv(x1,aa1)),prec);
  594.             t=gdiv(gaddsg(-1,t),gaddsg(1,t));
  595.           }
  596.         else t=gaddsg(2,ginv(gmul(u2,x1)));
  597.         tetpil=avma;return gerepile(av,tetpil,gcopy(t));
  598.       }
  599.     else err(zeller1);
  600. }  
  601.  
  602. GEN pointell(e,z,prec)
  603.      GEN e,z;
  604.      long prec;
  605. {
  606.   long av=avma,tetpil;
  607.   GEN p1,pii2,pii4,pii6,a,tau,q,u,y,yp,u1,u2,qn,qnu2,qnu1,qnu,qnu4,qnu3,p2,v;
  608.  
  609.   if((typ(e)!=17)||(lg(e)<20)) err(elliper1);
  610.   p1=mppi(prec);setexpo(p1,2);
  611.   pii2=cgetg(3,6);pii2[1]=zero;pii2[2]=(long)p1;
  612.   z=gdiv(z,e[15]);tau=gdiv(e[16],e[15]);
  613.   a=ground(gdiv(gimag(z),gimag(tau)));z=gsub(z,gmul(a,tau));
  614.   a=ground(greal(z));z=gsub(z,a);
  615.   if(gcmp0(z)||gexpo(gnorm(z))< -((prec-2)<<6)+10) {avma=av;v=cgetg(2,17);v[1]=zero;return v;}
  616.   q=gexp(gmul(pii2,tau),prec);u=gexp(gmul(pii2,z),prec);
  617.   y=gadd(gdivgs(gun,12),gdiv(u,u2=gsqr(u1=gsubsg(1,u))));
  618.   yp=gdiv(gaddsg(1,u),gmul(u1,u2));
  619.   qn=q;pii2=gdiv(pii2,e[15]);pii4=gmul(pii2,pii2);pii6=gmul(pii4,pii2);
  620.   do
  621.     {
  622.       p1=gsub(gmul(u,gadd(gdivsg(1,qnu2=gsqr(qnu1=gsubsg(1,qnu=gmul(qn,u)))),gdivsg(1,qnu4=gsqr(qnu3=gsub(qn,u))))),gmul2n(gdivsg(1,gsqr(gsubsg(1,qn))),1));
  623.       y=gadd(y,gmul(qn,p1));
  624.       p2=gadd(gdiv(gaddsg(1,qnu),gmul(qnu1,qnu2)),gdiv(gadd(qn,u),gmul(qnu3,qnu4)));
  625.       yp=gadd(yp,gmul(qn,p2));qn=gmul(q,qn);
  626.     }
  627.   while(gexpo(qn)>-((prec-2)<<5)-5);
  628.   yp=gmul(u,gmul(pii6,yp));y=gsub(gmul(pii4,y),gdivgs(e[6],12));
  629.   yp=gsub(yp,gadd(e[3],gmul(e[1],y)));
  630.   tetpil=avma;v=cgetg(3,17);v[1]=lcopy(y);v[2]=lmul2n(yp,-1);
  631.   return gerepile(av,tetpil,v);
  632. }
  633.  
  634.  
  635. GEN apell2(e,p)
  636.      GEN e,p; 
  637. {
  638. /* Calcul de a_p par les sommes de Jacobi */
  639.   GEN y,e71;
  640.   long av=avma,av2,i,l,s,e72,e6,e8;
  641.  
  642.   if((typ(e)!=17)||(lg(e)<14)) err(elliper1);
  643.   if(lgef(p)>3) err(apeller1);
  644.   s=0;l=p[2];e71=gmul2n(e[7],1);av2=avma;
  645.   if((unsigned long)l>=0x40000000) err(apeller1);
  646.   if((l<757)&&(l>2))
  647.     {
  648.       e6=itos(modis(e[6],l));e72=(itos(modis(e[7],l)))<<1;
  649.       e8=itos(modis(e[8],l));
  650.       for(i=0;i<l;i++) s=s+kross(e8+i*(e72+i*(e6+(i<<2))),l);
  651.     }
  652.   else
  653.     {
  654.       for(i=0;i<l;i++)
  655.         {
  656.           y=gadd(e[8],gmulsg(i,gadd(e71,gmulsg(i,gaddsg(i<<2,e[6])))));
  657.           s=s+krogs(y,l);avma=av2;
  658.         }
  659.     }
  660.   avma=av;return stoi(-s);
  661. }
  662.  
  663. GEN addsell(e,z1,z2)
  664.      GEN e,z1,z2;
  665.   
  666. {
  667.   GEN y,p1,p2,x1,x2,x3,y1,y2,y3,al;
  668.   long av=avma,tetpil;
  669.  
  670.   if(lg(z1)<3) return gcopy(z2);
  671.   if(lg(z2)<3) return gcopy(z1);
  672.   x1=(GEN)z1[1];x2=(GEN)z2[1];y1=(GEN)z1[2];y2=(GEN)z2[2];
  673.   if(gegal(x1,x2))
  674.     {
  675.       if(!gegal(y1,y2)) {y=cgetg(2,17);y[1]=zero;return y;}
  676.       else
  677.         {
  678.           p1=gadd(e,gmul(x1,gmulsg(3,x1)));
  679.           p2=gmul2n(y1,1);
  680.           if(gcmp0(p2)) {avma=av;y=cgetg(2,17);y[1]=zero;return y;}
  681.         }
  682.     }
  683.   else {p1=gsub(y2,y1);p2=gsub(x2,x1);}
  684.   al=gdiv(p1,p2);
  685.   x3=gsub(gmul(al,al),gadd(x1,x2));
  686.   y3=gneg(gadd(y1,gmul(al,gsub(x3,x1))));
  687.   tetpil=avma;y=cgetg(3,17);y[1]=lcopy(x3);y[2]=lcopy(y3);
  688.   return gerepile(av,tetpil,y);
  689. }
  690.  
  691. GEN doubsell(e,z1)
  692.      GEN e,z1;
  693. {
  694.   GEN x1,x3,y3,y,y1,p1,p2,al;
  695.   long av=avma,tetpil;
  696.  
  697.   if(lg(z1)<3) return gcopy(z1);
  698.   x1=(GEN)z1[1];y1=(GEN)z1[2];
  699.   p1=gadd(e,gmul(x1,gmulsg(3,x1)));p2=gmul2n(y1,1);
  700.   if(gcmp0(p2)) {avma=av;y=cgetg(2,17);y[1]=zero;return y;}
  701.   al=gdiv(p1,p2);
  702.   x3=gsub(gmul(al,al),gadd(x1,x1));
  703.   y3=gneg(gadd(y1,gmul(al,gsub(x3,x1))));
  704.   tetpil=avma;y=cgetg(3,17);y[1]=lcopy(x3);y[2]=lcopy(y3);
  705.   return gerepile(av,tetpil,y);
  706. }
  707.  
  708. GEN subsell(e,z1,z2)
  709.      GEN e,z1,z2;
  710.  
  711. {
  712.   GEN zp;
  713.   long av=avma,tetpil;
  714.   
  715.   if(lg(z2)<3) return gcopy(z1);
  716.   zp=cgetg(3,17);zp[1]=z2[1];
  717.   zp[2]=lneg(z2[2]);
  718.   tetpil=avma;return gerepile(av,tetpil,addsell(e,z1,zp));
  719. }
  720.  
  721.  
  722. GEN powsell(e,n,z)
  723.      GEN e,n,z;
  724.  
  725. {
  726.   GEN y,zp;
  727.   long s=signe(n),av=avma,i,j,tetpil;
  728.   unsigned long m;
  729.  
  730.   if(!s) {y=cgetg(2,17);y[1]=zero;return y;}
  731.   if(lg(z)<3) return gcopy(z);
  732.   if(s<0) 
  733.     {
  734.       n=gneg(n);zp=cgetg(3,17);zp[1]=z[1];
  735.       zp[2]=lneg(z[2]);
  736.     }
  737.   else zp=z;
  738.   if(gcmp1(n)) {tetpil=avma;return gerepile(av,tetpil,gcopy(zp));}
  739.   y=cgetg(2,17);y[1]=zero;
  740.   for (i=lgef(n)-1;i>2;i--)
  741.     {
  742.       for (m=n[i],j=0;j<32;j++,m>>=1)
  743.         {
  744.           if (m&1) y=addsell(e,y,zp);
  745.           zp=doubsell(e,zp);
  746.         }
  747.     }
  748.   for (m=n[2];m>1;m>>=1)
  749.     {
  750.       if (m&1) y=addsell(e,y,zp);
  751.       zp=doubsell(e,zp);
  752.     }
  753.   tetpil=avma;y=addsell(e,y,zp);
  754.   return gerepile(av,tetpil,y);
  755. }
  756.  
  757. GEN powssell(e,n,z)
  758.      GEN e,z;
  759.      long n;
  760.  
  761. {
  762.   GEN y,zp;
  763.   long av=avma,tetpil;
  764.   unsigned long m;
  765.  
  766.   if(!n) {y=cgetg(2,17);y[1]=zero;return y;}
  767.   if(lg(z)<3) return gcopy(z);
  768.   if(n<0) 
  769.     {
  770.       n= -n;zp=cgetg(3,17);zp[1]=z[1];
  771.       zp[2]=lneg(z[2]);
  772.     }
  773.   else zp=z;
  774.   if(n==1) {tetpil=avma;return gerepile(av,tetpil,gcopy(zp));}
  775.   y=cgetg(2,17);y[1]=zero;
  776.   for (m=n;m>1;m>>=1)
  777.     {
  778.       if (m&1) y=addsell(e,y,zp);
  779.       zp=doubsell(e,zp);
  780.     }
  781.   tetpil=avma;y=addsell(e,y,zp);
  782.   return gerepile(av,tetpil,y);
  783. }
  784.  
  785. #define HASHSP 255
  786.  
  787. GEN apell1(e,p)
  788.      GEN e,p;
  789. {
  790.   long av,av3,tetpil,k,k2,i,j,j1,lim,limite,succes,com,j2,s;
  791.   GEN tabla, tablb, count, index, hash;
  792.   GEN p1,p2,p3,q,h,hp,f,fh,fg,ftest;
  793.   GEN unmodp,pordmin,u,p1p,p2p,acon,bcon,xp,yp,c4,c6,cp4;
  794.   long flc,flcc,sucfin,fll,flf,x,pfinal;
  795.  
  796.   if((typ(e)!=17)||(lg(e)<14)) err(elliper1);
  797.   if(gcmpgs(p,20)<0) return apell2(e,p);
  798.   tabla = newbloc(1000000);
  799.   tablb = newbloc(1000000);
  800.   count = newbloc(256);
  801.   index = newbloc(257);
  802.   hash = newbloc(1000000);
  803.  
  804.   av=avma;limite=(av+bot)>>1;
  805.   unmodp=gmodulcp(gun,p);c4=gdivgs(gmul(unmodp,e[10]),-48);
  806.   c6=gdivgs(gmul(unmodp,e[11]),-864);
  807.   pordmin=gceil(gmul2n(gsqrt(p,4),2));p2p=gmul2n(p1p=gaddsg(1,p),1);
  808.   x=0;flcc=0;flc=kronecker(c6[2],p);u=c6;acon=gzero;bcon=gun;
  809.   sucfin=1;h=p1p;
  810.   while(sucfin)
  811.     {
  812.       while((flc==flcc)||(!flc))
  813.         {
  814.           x++;u=gadd(c6,gmulsg(x,gaddgs(c4,x*x)));
  815.           flc=kronecker(u[2],p);
  816.         }
  817.       flcc=flc;
  818.       s=itos(gceil(gsqrt(gdiv(pordmin,bcon),3)))>>1;
  819.       cp4=gmul(c4,yp=gsqr(u));
  820.       xp=gmulsg(x,u);f=cgetg(3,17);f[1]=(long)xp;f[2]=(long)yp;
  821.       fh=powsell(cp4,h,f);
  822.       if (bcon != gun) f=powsell(cp4,bcon,f); /* sic */
  823.       p1=fh;flf=1;
  824.       for(i=0;i<=HASHSP;i++) count[i]=0;
  825.       for(i=0;(i<=s-1)&&flf;i++)
  826.         {
  827.           if(lg(p1)==3)
  828.             {
  829.              
  830. p2=(GEN)((GEN)p1[1])[2];tabla[i]=p2[lgef(p2)-1];j=tabla[i]&HASHSP;count[j]++;
  831. /*       printf("%d ",j);fflush(stdout); */
  832.               p2=(GEN)((GEN)p1[2])[2];tablb[i]=p2[lgef(p2)-1];
  833.               p1=addsell(cp4,p1,f);
  834.             }
  835.           else flf=0;
  836.         }
  837. /*      printf("\nsj:\n"); */
  838.       if(flf)
  839.         {
  840.           fg=powssell(cp4,s,f);ftest=fg;
  841.           index[0]=0;for(i=0;i<=HASHSP-1;i++) index[i+1]=index[i]+count[i];
  842.           for(i=0;i<=s-1;i++) hash[index[tabla[i]&HASHSP]++]=i;
  843.           index[0]=0;for(i=0;i<=HASHSP;i++) index[i+1]=index[i]+count[i];
  844.           succes=0;com=1;av3=avma;pfinal=p[lgef(p)-1];
  845.           while(!succes)
  846.             {
  847.               p1=(GEN)((GEN)ftest[1])[2];k=p1[lgef(p1)-1];j=k&HASHSP;
  848. /*      printf("%d ",j);fflush(stdout); */
  849.               fll=(lg(ftest)==3);
  850.               for(j1=index[j];(j1<index[j+1])&&(!succes);j1++)
  851.                 {
  852.                   j2=hash[j1];
  853.                   if((tabla[j2]==k)&&fll)
  854.                     {
  855.                       p2=(GEN)((GEN)ftest[2])[2];k2=p2[lgef(p2)-1];
  856.                       if((tablb[j2]==k2)||(tablb[j2]==pfinal-k2))
  857.                         {
  858.                           p1=addsell(cp4,powssell(cp4,j2,f),fh);
  859.                           succes=gegal(p1[1],ftest[1]);
  860.                         }
  861.                     }
  862.                 }
  863.               if(!succes)
  864.                 {
  865.                   com++;
  866.                   if(avma>=limite) ftest=addsell(cp4,ftest,fg);
  867.                   else 
  868.                     {
  869.                      
  870. tetpil=avma;ftest=gerepile(av3,tetpil,addsell(cp4,ftest,fg));
  871.                     }
  872.                   if (lg(ftest)<3) err(apeller2);
  873.                 }
  874.             }
  875.           h=addii(h,mulsi(j2,bcon));p2=mulsi(s,mulsi(com,bcon));
  876.          
  877. h=(!cmpii(((GEN)p1[2])[2],((GEN)ftest[2])[2]))?subii(h,p2):addii(h,p2);
  878.         }
  879.       else h=addii(h,mulsi(i-1,bcon));
  880.       p2=factor(h);
  881.       p1=(GEN)p2[1];
  882.       p2=(GEN)p2[2];
  883.       for(i=1;i<lg(p1);i++)
  884.         {
  885.           p3=divii(h,p1[i]);fh=powsell(cp4,p3,f);lim=itos(p2[i]);
  886.           for(j=1;(j<=lim)&&(lg(fh)<3);j++)
  887.             {
  888.               h=p3;if(j<lim) {p3=divii(h,p1[i]);fh=powsell(cp4,p3,f);}
  889.             }
  890.         }
  891.       p1=gmodulcp(acon,bcon);p2=gmodulcp(gzero,h);
  892.       p1=chinois(p1,p2);acon=(GEN)p1[2];bcon=(GEN)p1[1];
  893.       if(gcmp(bcon,pordmin)>=0)
  894.         {
  895.           q=ground(gdiv(gsub(p1p,acon),bcon));sucfin=0;
  896.           hp=addii(mulii(q,bcon),acon);tetpil=avma;
  897.         }
  898.       else
  899.         {
  900.           acon=modii(subii(p2p,acon),bcon);
  901.           p1=subii(acon,bcon);if(signe(addii(acon,p1))>0) acon=p1;
  902.           q=ground(gdiv(gsub(p1p,acon),bcon));
  903.           h=addii(mulii(q,bcon),acon);
  904.         }
  905.     }
  906.   killbloc(tabla); killbloc(tablb); killbloc(count);
  907.   killbloc(index); killbloc(hash);
  908.   return
  909. (flc==1)?gerepile(av,tetpil,gsub(p1p,hp)):gerepile(av,tetpil,gsub(hp,p1p));
  910. }
  911.  
  912. typedef struct {int isnull; long x; long y;} sellpt;
  913.  
  914. long addmod(a, b, p)
  915.      long a, b, p;
  916. {
  917.   long res = a + b;
  918.   return (res >= p) ? res - p : res;
  919. }
  920.  
  921. long submod(a, b, p)
  922.      long a, b, p;
  923. {
  924.   long res = a - b;
  925.   return (res >= 0) ? res : res + p;
  926. }
  927.  
  928. extern long mulmodll();
  929. #define mulmod(a,b,c) mulmodll(a,b,c)
  930.  
  931. static long divmod(a, b, p)
  932.      long a, b, p;
  933. {
  934.   long v1 = 0, v2 = 1, v3, r, oldp = p;
  935.  
  936.   while (b > 1)
  937.     {
  938.       v3 = v1 - (p / b) * v2; v1 = v2; v2 = v3;
  939.       r = p % b; p = b; b = r;
  940.     }
  941.  
  942.   if (v2 < 0) v2 += oldp;
  943.   return mulmod(a, v2, oldp);
  944. }
  945.   
  946. void addsell1(e, p, p1, p2, p3)
  947.      long e, p;
  948.      sellpt *p1, *p2, *p3;
  949. {
  950.   long num, den, lambda;
  951.   if (p1->isnull) {*p3 = *p2; return;}
  952.   if (p2->isnull) {*p3 = *p1; return;}
  953.   if (p1->x == p2->x)
  954.     if ((p1->y)&&(p1->y == p2->y))
  955.       {
  956.         num = addmod(e, mulmod(3, mulmod(p1->x, p1->x, p), p), p);
  957.         den = addmod(p1->y, p1->y, p);
  958.       }
  959.     else
  960.       {
  961.         p3->isnull = 1;
  962.         return;
  963.       }
  964.   else
  965.     {
  966.       num = submod(p1->y, p2->y, p);
  967.       den = submod(p1->x, p2->x, p);
  968.     }
  969.   lambda = divmod(num, den, p);
  970.   p3->x = submod(mulmod(lambda, lambda, p), addmod(p1->x, p2->x, p), p);
  971.   p3->y = submod(mulmod(lambda, submod(p2->x, p3->x, p), p), p2->y, p);
  972.   p3->isnull = 0;
  973. }
  974.  
  975. void powssell1(e, p, n, p1, p2)
  976.      long e, p, n;
  977.      sellpt *p1, *p2;
  978. {
  979.   sellpt p4, p3;
  980.   p3 = *p1;
  981.   if (n < 0) {n = -n; if (p3.y) p3.y = p - p3.y;}
  982.   p2->isnull = 1;
  983.   for(;;)
  984.     {
  985.       if (n&1) addsell1(e, p, p2, &p3, p2);
  986.       n>>=1;
  987.       if (!n) return;
  988.       addsell1(e, p, &p3, &p3, &p4);
  989.       p3 = p4;
  990.     }
  991. }
  992.  
  993. typedef struct {long x; long y; long i;} multiple;
  994.  
  995. int compare_multiples(a, b)
  996.      multiple *a, *b;
  997. {
  998.   return a->x - b->x;
  999. }
  1000.  
  1001. GEN apell(e,pl)
  1002.      GEN e,pl;
  1003. {
  1004.   long av = avma,i,j,com,s;
  1005.   GEN p1,p2,unmodp;
  1006.   sellpt f,fh,fg,ftest,f2;
  1007.   long pordmin,u,p1p,p2p,acon,bcon,xp,yp,c4,c6,cp4;
  1008.   long flb,flc,flcc,x, q, h, p3, p, l, r, m;
  1009.   multiple *table;
  1010.  
  1011.   if((typ(e)!=17)||(lg(e)<14)) err(elliper1);
  1012.   if(divise(e[12],pl)) return stoi(((((pl[lgef(pl)-1]&3)+1)&2)-1)*kronecker(e[11],pl));
  1013.   if (gcmpgs(pl, 0x3fffffff) > 0) return apell1(e, pl);
  1014.   p = itos(pl);
  1015.   if (p < 100) return apell2(e, pl);
  1016.   unmodp = gmodulcp(gun, pl);
  1017.   c4 = itos(gdivgs(gmul(unmodp, e[10]), -48)[2]);
  1018.   c6 = itos(gdivgs(gmul(unmodp, e[11]), -864)[2]);
  1019.   pordmin = 1 + 4 * sqrt((float)p);
  1020.   p1p = p + 1; p2p = p1p << 1;
  1021.   x = 0; flcc = 0; flc = kross(c6, p); u = c6; acon = 0; bcon = 1;
  1022.   h = p1p;
  1023.   for(;;)
  1024.     {
  1025.       while((flc == flcc) || (!flc))
  1026.         {
  1027.           x++;
  1028.           u = addmod(c6, mulmod(x, c4 + mulmod(x, x, p), p), p);
  1029.           flc = kross(u,p);
  1030.         }
  1031.       flcc = flc;
  1032.       s = sqrt(((float)pordmin)/bcon) / 2;
  1033.       if (!s) s = 1;
  1034.       if(bcon==1) table=(multiple *)malloc((s+1)*sizeof(multiple));
  1035.       yp = mulmod(u, u, p);
  1036.       cp4 = mulmod(c4, yp, p);
  1037.       xp = mulmod(x, u, p);
  1038.       f.isnull = 0; f.x = xp; f.y = yp;
  1039.       powssell1(cp4, p, h, &f, &fh);
  1040.       if (bcon > 1) powssell1(cp4, p, bcon, &f, &f2);else f2=f;
  1041.       for(i=0; i < s; i++)
  1042.         if (fh.isnull)
  1043.           {
  1044.             h += bcon * i;
  1045.             goto trouve;
  1046.           }
  1047.         else
  1048.           {
  1049.             table[i].x = fh.x;
  1050.             table[i].y = fh.y;
  1051.             table[i].i = i;
  1052.             addsell1(cp4, p, &fh, &f2, &fh);
  1053.           }
  1054.       qsort(table, s, sizeof(multiple), compare_multiples);
  1055.       powssell1(cp4, p, s, &f2, &fg); ftest = fg;
  1056.       for(com = 1;; com++)
  1057.         {
  1058.           if (ftest.isnull) err(apeller3);
  1059.           l = 0; r = s;
  1060.           while (l < r)
  1061.             {
  1062.               m = (l + r) >> 1;
  1063.               if (table[m].x < ftest.x) l = m + 1; else r = m;
  1064.             }
  1065.           if ((r < s) && (table[r].x == ftest.x)) break;
  1066.           addsell1(cp4, p, &ftest, &fg, &ftest);
  1067.         }
  1068.       h += table[r].i * bcon;
  1069.       if (table[r].y == ftest.y) h -= s * com * bcon; else h += s * com *
  1070. bcon;
  1071.  
  1072.     trouve:
  1073.       p2=factor(stoi(h));
  1074.       p1=(GEN)p2[1];
  1075.       p2=(GEN)p2[2];
  1076.       for(i=1; i < lg(p1); i++)
  1077.         for(j = ((GEN)p2[i])[2]; j; j--)
  1078.         {
  1079.           p3 = h / ((GEN)p1[i])[2];
  1080.           powssell1(cp4, p, p3, &f, &fh);
  1081.           if (!fh.isnull) break;
  1082.           h = p3;
  1083.         }
  1084.       flb=0;
  1085.       if (bcon > 1)
  1086.         {
  1087.           p1 = gmodulcp(stoi(acon), stoi(bcon)); p2=gmodulcp(gzero, stoi(h));
  1088.           p1=chinois(p1,p2);acon=itos((GEN)p1[2]);bcon=((GEN)p1[1])[2];
  1089.       if((bcon<0)||(lgef(p1[1])>3)) flb=1;
  1090.         }
  1091.       else
  1092.         bcon = h;
  1093.       if(flb||(bcon >= pordmin))
  1094.         {
  1095.           if(flb) h=acon;
  1096.       else
  1097.         {
  1098.           q = ((unsigned long)(p2p + bcon - (acon << 1))) / (bcon << 1);
  1099.           h = acon + q * bcon;
  1100.         }
  1101.           break;
  1102.         }
  1103.       else
  1104.         {
  1105.           acon = (p2p - acon) % bcon;
  1106.           if ((acon << 1) > bcon) acon -= bcon;
  1107.           q = ((unsigned long)(p2p + bcon - (acon << 1))) / (bcon << 1);
  1108.           h = acon + q * bcon;
  1109.         }
  1110.     }
  1111.   avma = av;free(table);
  1112.   return stoi((flc == 1) ?  p1p - h : h - p1p);
  1113. }
  1114.  
  1115. GEN anell(e, n)
  1116.      GEN e;
  1117.      long n;
  1118. {
  1119.   long p, pk, i, m, av, tetpil, pl[3];
  1120.   GEN p1, p2, ap, an;
  1121.   
  1122.   if((typ(e)!=17)||(lg(e)<14)) err(elliper1);
  1123.   if (n <= 0) return cgetg(1, 17);
  1124.   an = cgetg(n+1, 17);
  1125.   an[1] = un;
  1126.   for(i=2; i <= n; i++) an[i] = 0;
  1127.   pl[0] = 0x1000003; pl[1] = 0x1010003;
  1128.   for (p = 2; p <= n; p++) if (!an[p])
  1129.     {
  1130.       pl[2] = p;
  1131.       if (divise(e[12], pl)) /* mauvaise reduction */
  1132.         switch (((((p & 3) + 1) & 2) -1) * krogs(e[11], p)) /* renvoie (-c6 /
  1133. p) */
  1134.           {
  1135.           case -1:  /* non deployee */
  1136.             for(m = p; m <= n; m += p) if (an[m / p]) an[m] = lneg(an[m /
  1137. p]);
  1138.             continue;
  1139.           case 0:   /* additive */
  1140.             for(m = p; m <= n; m += p) an[m] = zero;
  1141.             continue;
  1142.           case 1:   /* deployee */
  1143.             for(m = p; m <= n; m += p) if (an[m / p]) an[m] = lcopy(an[m /
  1144. p]);
  1145.           }
  1146.       else /* bonne reduction */
  1147.         for(pk = p; pk <= n; pk *= p)
  1148.           {
  1149.             if (pk == p)
  1150.               an[pk] = (long)(ap = apell(e, pl));
  1151.             else
  1152.               {
  1153.                 av = avma;
  1154.                 p1 = mulii(ap, an[pk / p]);
  1155.                 p2 = mulsi(p, an[pk / p / p]);
  1156.                 tetpil = avma;
  1157.                 an[pk] = lpile(av, tetpil, subii(p1, p2));
  1158.               }
  1159.             for(m = n / pk; m > 1; m--)
  1160.               if (an[m] && (m % p)) an[m * pk] = lmulii(an[m], an[pk]);
  1161.           }
  1162.     }
  1163.   return an;
  1164. }
  1165.  
  1166. GEN akell(e, n)
  1167.      GEN e,n;
  1168. {
  1169.   long i,j,ex,av=avma,tetpil;
  1170.   GEN p1,p2,ap,u,v,w,fac,y,pl;
  1171.   
  1172.   if((typ(e)!=17)||(lg(e)<14)) err(elliper1);
  1173.   if(typ(n)!=1) err(akeller1);
  1174.   if(signe(n)<= 0) return gzero;
  1175.   y=gun;if(gcmp1(n)) return y;
  1176.   fac=auxdecomp(n,1);p1=(GEN)fac[1];p2=(GEN)fac[2];
  1177.   for(i=1;(i<lg(p1))&&(!gcmp0(y));i++)
  1178.     {
  1179.       pl=(GEN)p1[i];
  1180.       if (divise(e[12], pl)) /* mauvaise reduction */
  1181.         {
  1182.       j=((((pl[lgef(pl)-1]&3)+1)&2)-1)*kronecker(e[11],pl);
  1183.       if((j<0)&&(mpodd(p2[i]))) {tetpil=avma;y=gneg(y);}
  1184.       if(!j) {avma=av;y=gzero;}
  1185.     }
  1186.       else /* bonne reduction */
  1187.     {
  1188.       ap=apell(e,pl);ex=itos(p2[i]);
  1189.       u=ap;v=gun;
  1190.       for(j=2;j<=ex;j++)
  1191.         {
  1192.           w=gsub(gmul(ap,u),gmul(pl,v));
  1193.           v=u;u=w;
  1194.         }
  1195.       tetpil=avma;y=gmul(u,y);
  1196.     }
  1197.     }
  1198.   return (av==avma)?y:gerepile(av,tetpil,y);
  1199. }
  1200.  
  1201.  
  1202. GEN hell3(e,a,prec)
  1203.      GEN e,a;
  1204.      long prec;
  1205. {
  1206.   long av=avma,tetpil;
  1207.   GEN p1,p2,z,q,w,piisurw,pitemp;
  1208.  
  1209.   if((typ(e)!=17)||(lg(e)<20)) err(elliper1);
  1210.   if(!oncurve(e,a)) err(heller1);
  1211.   if(lg(a)<3) return gzero;
  1212.   constpi(prec);
  1213.   z=zell(e,a,prec);pitemp=mppi(prec);piisurw=gdiv(gmul(gi,pitemp),w=(GEN)e[15]);
  1214.   q=gexp(gmul(e[16],piisurw),prec);
  1215.   p1=gdiv(theta(q,gmul(gdiv(z,w),pitemp),prec),thetanullk(q,1,prec));
  1216.   p1=gmul(p1,gdiv(w,pitemp));
  1217.   p1=glog(gmul(p1,gdiv(denom(a[1]),denom(a[2]))),prec);
  1218.   p2=gimag(z);
  1219.   if((gcmp0(p2))||((gexpo(p2)-gexpo(gimag(e[16])))<-2))
  1220.     p1=gneg(p1);
  1221.   else {p2=gmul(gmul2n(pitemp,-2),gimag(gdiv(e[16],e[15])));p1=gsub(p2,p1);}
  1222.   tetpil=avma;return gerepile(av,tetpil,greal(p1));
  1223. }
  1224.  
  1225. GEN hell(e,a,prec)
  1226.      GEN e,a;
  1227.      long prec;
  1228. {
  1229.   long av=avma,tetpil,n;
  1230.   GEN p1,p2,y,z,q,psi2,pi2surw,pi2isurw,qn,ps;
  1231.  
  1232.   if((typ(e)!=17)||(lg(e)<20)) err(elliper1);
  1233.   pi2surw=gdiv(gmul2n(mppi(prec),1),(GEN)e[15]);
  1234.   pi2isurw=cgetg(3,6);pi2isurw[1]=zero;pi2isurw[2]=(long)pi2surw;
  1235.   z=gmul(greal(zell(e,a,prec)),pi2surw);
  1236.   q=greal(gexp(gmul(e[16],pi2isurw),prec));
  1237.   psi2=gadd(e[3],gadd(gmul(e[1],a[1]),gmul2n(a[2],1)));
  1238.   y=gsin(z,prec);n=0;qn=gun;ps=gneg(q);
  1239.   do
  1240.     {
  1241.       n++;p1=gsin(gmulsg(n+n+1,z),prec);qn=gmul(qn,ps);
  1242.       ps=gmul(ps,q);p1=gmul(p1,qn);
  1243.       y=gadd(y,p1);
  1244.     }
  1245.   while(gexpo(qn)>= -((prec-2)<<5));
  1246.   p1=gmul(gsqr(gdiv(gmul2n(y,1),psi2)),pi2surw);
  1247.   p2=gsqr(gsqr(gdiv(p1,gsqr(gsqr(denom(a[1]))))));
  1248.   p1=gdiv(gmul(p2,q),e[12]);
  1249.   p1=gmul2n(glog(gabs(p1,prec),prec),-5);
  1250.   tetpil=avma;return gerepile(av,tetpil,gneg(p1));
  1251. }
  1252.  
  1253. GEN ghell(e,a,prec)
  1254.      GEN e,a;
  1255.      long prec;
  1256. {
  1257.   long av=avma,av2,tetpil,lx,i,n,n2,grandn;
  1258.   GEN p,p1,p2,x,y,z,phi2,psi2,psi3,logdep;
  1259.  
  1260. /* On suppose que la courbe est a coefficients entiers donne par un
  1261. modele minimal */
  1262.  
  1263.   if((typ(e)!=17)||(lg(e)<20)) err(elliper1);
  1264.   if(!oncurve(e,a)) err(heller1);
  1265.   if(lg(a)<3) return gzero;
  1266.   x=(GEN)a[1];y=(GEN)a[2];
  1267.   psi2=numer(gadd(e[3],gadd(gmul(e[1],x),gmul2n(y,1))));
  1268.   p2=gadd(gmulsg(3,e[7]),gmul(x,gadd(e[6],gmulsg(3,x))));
  1269.   psi3=numer(gadd(e[9],gmul(x,gadd(gmulsg(3,e[8]),gmul(x,p2)))));
  1270.   if((!signe(psi2))||(!signe(psi3))) {avma=av;return gzero;}
  1271.  
  1272. phi2=numer(gsub(gadd(e[4],gmul(x,gadd(shifti(e[2],1),gmulsg(3,x)))),gmul(e[1],
  1273. y)));
  1274.   p1=(GEN)factor(mppgcd(psi2,phi2))[1];lx=lg(p1);
  1275.   tetpil=avma;z=hell(e,a,prec);av2=avma;
  1276.   for(i=1;i<lx;i++)
  1277.     {
  1278.       p=(GEN)p1[i];n2=ggval(psi2,p);logdep=gneg(glog(p,prec));
  1279.       if(signe(resii(e[10],p)))
  1280.         {
  1281.           grandn=ggval(e[12],p);n=n2<<1;
  1282.       if(grandn)
  1283.         {
  1284.           if(n>grandn) n=grandn;
  1285.           p2=divrs(mulsr(n*(grandn+grandn-n),logdep),(grandn<<3));
  1286.           tetpil=avma;z=gadd(z,p2);
  1287.         }
  1288.       else avma=av2;
  1289.         }
  1290.       else
  1291.         {
  1292.           n=ggval(psi3,p);
  1293.           if(n>=3*n2) p2=gdivgs(mulsr(n2,logdep),3);
  1294.           else p2=gmul2n(mulsr(n,logdep),-3);
  1295.           tetpil=avma;z=gadd(z,p2);
  1296.         }
  1297.     }
  1298.   return gerepile(av,tetpil,z);
  1299. }
  1300.  
  1301. GEN ghell2(e,a,prec)
  1302.      GEN e,a;
  1303.      long prec;
  1304. {
  1305.   long av=avma,tetpil,av2,lx,i,n,n2,grandn;
  1306.   GEN p,p1,p2,x,y,z,phi2,psi2,psi3,logdep;
  1307.  
  1308. /* On suppose que la courbe est a coefficients entiers donne par un
  1309. modele minimal */
  1310.  
  1311.   if((typ(e)!=17)||(lg(e)<20)) err(elliper1);
  1312.   if(!oncurve(e,a)) err(heller1);
  1313.   if(lg(a)<3) return gzero;
  1314.   x=(GEN)a[1];y=(GEN)a[2];
  1315.   psi2=numer(gadd(e[3],gadd(gmul(e[1],x),gmul2n(y,1))));
  1316.   if(!signe(psi2)) return gzero;
  1317.  
  1318. phi2=numer(gsub(gadd(e[4],gmul(x,gadd(shifti(e[2],1),gmulsg(3,x)))),gmul(e[1],
  1319. y)));
  1320.   p1=(GEN)factor(mppgcd(psi2,phi2))[1];lx=lg(p1);
  1321.   tetpil=avma;z=hell2(e,a,prec);av2=avma;
  1322.   if(lx==1) return gerepile(av,tetpil,z);
  1323.   for(i=1;i<lx;i++)
  1324.     {
  1325.       p=(GEN)p1[i];n2=ggval(psi2,p);logdep=gneg(glog(p,prec));
  1326.       if(signe(resii(e[10],p)))
  1327.         {
  1328.           grandn=ggval(e[12],p);n=n2<<1;
  1329.       if(grandn)
  1330.         {
  1331.           if(n>grandn) n=grandn;
  1332.           p2=divrs(mulsr(n*(grandn+grandn-n),logdep),(grandn<<3));
  1333.           tetpil=avma;z=gadd(z,p2);
  1334.         }
  1335.       else avma=av2;
  1336.         }
  1337.       else
  1338.         {
  1339.           p2=gadd(gmulsg(3,e[7]),gmul(x,gadd(e[6],gmulsg(3,x))));
  1340.           psi3=numer(gadd(e[9],gmul(x,gadd(gmulsg(3,e[8]),gmul(x,p2)))));
  1341.           n=ggval(psi3,p);
  1342.           if(n>=3*n2) p2=gdivgs(mulsr(n2,logdep),3);
  1343.           else p2=gmul2n(mulsr(n,logdep),-3);
  1344.           tetpil=avma;z=gadd(z,p2);
  1345.         }
  1346.     }
  1347.   return gerepile(av,tetpil,z);
  1348. }
  1349.  
  1350. GEN ghell3(e,a,prec)
  1351.      GEN e,a;
  1352.      long prec;
  1353. {
  1354.   long av=avma,tetpil,av2,lx,i,n,n2,grandn;
  1355.   GEN p,p1,p2,x,y,z,phi2,psi2,psi3,logdep;
  1356.  
  1357. /* On suppose que la courbe est a coefficients entiers donne par un
  1358. modele minimal */
  1359.  
  1360.   if((typ(e)!=17)||(lg(e)<20)) err(elliper1);
  1361.   if(!oncurve(e,a)) err(heller1);
  1362.   if(lg(a)<3) return gzero;
  1363.   x=(GEN)a[1];y=(GEN)a[2];
  1364.   psi2=numer(gadd(e[3],gadd(gmul(e[1],x),gmul2n(y,1))));
  1365.   if(!signe(psi2)) return gzero;
  1366.  
  1367. phi2=numer(gsub(gadd(e[4],gmul(x,gadd(shifti(e[2],1),gmulsg(3,x)))),gmul(e[1],
  1368. y)));
  1369.   p1=(GEN)factor(mppgcd(psi2,phi2))[1];lx=lg(p1);
  1370.   tetpil=avma;z=hell3(e,a,prec);av2=avma;
  1371.   if(lx==1) return gerepile(av,tetpil,z);
  1372.   for(i=1;i<lx;i++)
  1373.     {
  1374.       p=(GEN)p1[i];n2=ggval(psi2,p);logdep=gneg(glog(p,prec));
  1375.       if(signe(resii(e[10],p)))
  1376.         {
  1377.           grandn=ggval(e[12],p);n=n2<<1;
  1378.       if(grandn)
  1379.         {
  1380.           if(n>grandn) n=grandn;
  1381.           p2=divrs(mulsr(n*(grandn+grandn-n),logdep),(grandn<<3));
  1382.           tetpil=avma;z=gadd(z,p2);
  1383.         }
  1384.       else avma=av2;
  1385.         }
  1386.       else
  1387.         {
  1388.           p2=gadd(gmulsg(3,e[7]),gmul(x,gadd(e[6],gmulsg(3,x))));
  1389.           psi3=numer(gadd(e[9],gmul(x,gadd(gmulsg(3,e[8]),gmul(x,p2)))));
  1390.           n=ggval(psi3,p);
  1391.           if(n>=3*n2) p2=gdivgs(mulsr(n2,logdep),3);
  1392.           else p2=gmul2n(mulsr(n,logdep),-3);
  1393.           tetpil=avma;z=gadd(z,p2);
  1394.         }
  1395.     }
  1396.   return gerepile(av,tetpil,z);
  1397. }
  1398.  
  1399. GEN lseriesell(e,s,N,A,prec)
  1400.      GEN e,s,N,A;
  1401.      long prec;
  1402. {
  1403.   long av=avma,tetpil,l,n,eps;
  1404.   GEN z,p1,p2,cg,cg1,v,cga,cgb,s2,ns,gs;
  1405.  
  1406.   if((typ(e)!=17)||(lg(e)<14)||(typ(N)!=1)||(!signe(N))) err(elliper1);
  1407.   if(gsigne(A)<=0) err(lserieser1);
  1408.   if(gcmpgs(A,1)<0) A=ginv(A);
  1409.   eps=signe(N);if(eps<0) N=gneg(N);
  1410.   cg1=mppi(prec);setexpo(cg1,2);cg=divrr(cg1,gsqrt(N,prec));
  1411.   cga=gmul(cg,A);cgb=gdiv(cg,A);
  1412.   l=(C2*(prec-2)+fabs(gtodouble(s)-1.)*log(rtodbl(cga)))/rtodbl(cgb)+1;
  1413.   v=anell(e,l);
  1414.   s2=gsubsg(2,s);ns=gpui(cg,gsubgs(gmul2n(s,1),2),prec);
  1415.   z=gzero;
  1416.   if(typ(s)==1)
  1417.     {
  1418.       if(signe(s)>0) gs=mpfactr(itos(s)-1,prec);
  1419.       else {avma=av;return gzero;}
  1420.     }
  1421.   else gs=ggamma(s,prec);/* gs2=ggamma(s2,prec); */
  1422.   for(n=1;n<=l;n++)
  1423.     {
  1424.       p1=gdiv(incgam4(s,gmulsg(n,cga),gs,prec),gpui(stoi(n),s,prec));
  1425.       p2=gdiv(gmul(ns,incgam(s2,gmulsg(n,cgb),prec)),gpui(stoi(n),s2,prec));
  1426.       if(eps<0) p2=gneg(p2);
  1427.       z=gadd(z,gmul(gadd(p1,p2),v[n]));
  1428.     }
  1429.   tetpil=avma;return gerepile(av,tetpil,gdiv(z,gs));
  1430. }
  1431.     
  1432. /********************************************************************/
  1433. /********************************************************************/
  1434. /**                                                                **/
  1435. /**                     Algorithme de Tate                         **/
  1436. /**                       (cf Anvers IV)                           **/
  1437. /**             Type de Kodaira, modele minimal global             **/
  1438. /**                                                                **/
  1439. /********************************************************************/
  1440. /********************************************************************/
  1441.  
  1442. /*
  1443.   Etant donnes une courbe elliptique sous forme longue e, dont les coefficients
  1444.   sont entiers, et un nombre premier p1, renvoie le type de la fibre en p du
  1445.   modele de Neron de la courbe, ainsi que les changements de variables 
  1446.   necessaires, sous la forme d'un triplet [f, kod, v].
  1447.   
  1448.   L'entier f est l'exposant du conducteur.
  1449.   
  1450.   l'entier kod, est le type de kodaira. les types II, III et IV sont codes 2,
  1451.   3, et 4 respectivement. Les types II*, III* et IV* donnent -2, -3 et -4. Le
  1452.   type I0* donne -1, les types Inu et Inu* donnent 4+nu et -4-nu. Enfin, le
  1453.   type I0 donne 1.
  1454.   
  1455.   v est un quadruplet [u, r, s, t] qui permet de passer a un modele minimal.
  1456.   En general, on ne s'interessera a ce vecteur que si u <> 1.
  1457.   
  1458.   L'algorithme est bien sur celui de Tate dans Anvers IV. Compte tenu des
  1459.   remarques du bas de la page 46, l'algorithme long n'est utilise que pour 
  1460.   p = 2 ou p = 3.
  1461.   
  1462.   Remarque. L'algorithme donne aussi le nombre c (voir Tate), qui n'est pas
  1463.   calcule ici. Peut-etre faudrait-il le faire ?
  1464.  
  1465.  */
  1466.  
  1467. static void cumule(), cumule1();
  1468.  
  1469. GEN localreduction(e, p1)
  1470.     GEN e, p1;
  1471. {
  1472.   long av = avma, tetpil, k, p, f, kod, nu, nuj, nudelta;
  1473.   int a21, a42, a63, a32, a64, theroot;
  1474.   GEN pk, p2k, pk1, p4, p6, a2prime, a3prime;
  1475.   GEN p2, p3, r = gzero, s = gzero, t = gzero, v, result;
  1476.   
  1477.   if ((typ(e) != 17) || (lg(e) < 14)) err(elliper1);
  1478.   
  1479.   v = cgetg(5,17); v[1] = un; v[2] = v[3] = v[4] = zero;
  1480.   nudelta = ggval(e[12], p1);
  1481.   if (gcmpgs(p1, 3) > 0)       /* p different de 2 ou 3 */
  1482.     {
  1483.       nuj = gcmp0(e[13]) ? 0 : - ggval(e[13], p1);
  1484.       k = (nuj > 0 ? nudelta - nuj : nudelta) / 12;
  1485.       if (k > 0) /* modele non minimal */
  1486.         {
  1487.           pk = gpuigs(p1, k);
  1488.           if (mpodd(e[1]))
  1489.             s = shifti(subii(pk, e[1]), -1);
  1490.           else
  1491.             s = negi(shifti(e[1], -1));
  1492.           p2k = mulii(pk, pk);
  1493.           a2prime = subii(e[2], mulii(s, addii(e[1], s)));
  1494.           switch(itos(modis(a2prime, 3)))
  1495.             {
  1496.             case 0: r = negi(divis(a2prime, 3)); break;
  1497.             case 1: r = divis(subii(p2k, a2prime), 3); break;
  1498.             case 2: r = negi(divis(addii(a2prime, p2k), 3)); break;
  1499.             default: err(tater1);
  1500.             }
  1501.           a3prime = addii(e[3], mulii(r, e[1]));
  1502.           if (mpodd(a3prime))
  1503.             t = shifti(subii(mulii(pk, p2k), a3prime), -1);
  1504.           else
  1505.             t = negi(shifti(a3prime, -1));
  1506.           v[1] = (long)pk; v[2] = (long)r; v[3] = (long)s; v[4] = (long)t;
  1507.           nudelta -= 12 * k;
  1508.         }
  1509.       if (nuj > 0) switch(nudelta - nuj)
  1510.         {
  1511.         case 0: f = 1; kod = 4 + nuj; break;               /* Inu */
  1512.         case 6: f = 2; kod = - 4 - nuj; break;             /* Inu* */
  1513.         default: err(tater1);
  1514.         }
  1515.       else switch(nudelta)
  1516.         {
  1517.         case 0: f = 0; kod = 1; break;                    /* I0, regulier */
  1518.         case 2: f = 2; kod = 2; break;                    /* II   */
  1519.         case 3: f = 2; kod = 3; break;                    /* III  */
  1520.         case 4: f = 2; kod = 4; break;                    /* IV   */
  1521.         case 6: f = 2; kod = -1; break;                   /* I0*  */
  1522.         case 8: f = 2; kod = -4; break;                   /* IV*  */
  1523.         case 9: f = 2; kod = -3; break;                   /* III* */
  1524.         case 10: f = 2; kod = -2; break;                  /* II*  */
  1525.         default: err(tater1);
  1526.         }
  1527.     }
  1528.   else for(;;)
  1529.     {
  1530.       /* ici, p = 2 ou p = 3 */
  1531.       if (!nudelta) {f = 0; kod = 1; goto exit;}                             
  1532.  /* I0   */
  1533.       if (!divise(e[6], p1)) {f =  1; kod = 4 + nudelta; goto exit;}         
  1534.  /* Inu  */
  1535.       p = itos(p1);
  1536.       if (p == 2)
  1537.         {
  1538.           r = modis(e[4], 2);
  1539.           s = modis(addii(r, e[2]), 2);
  1540.           if (signe(r)) t = modis(addii(addii(e[4], e[5]), s), 2);
  1541.           else t = modis(e[5], 2);
  1542.         }
  1543.       else /* p == 3 */
  1544.         {
  1545.           r = negi(modis(e[8], 3));
  1546.           s = modis(e[1], 3);
  1547.           t = modis(addii(e[3], mulii(e[1], r)), 3);
  1548.         }
  1549.       cumule(&v, &e, gun, r, s, t);       /* p | a1, a2, a3, a4 et a6 */
  1550.       p2 = stoi(p*p);
  1551.       if(!divise(e[5], p2)) {f = nudelta; kod = 2; goto exit;}               
  1552.   /* II   */
  1553.       p3 = stoi(p*p*p);
  1554.       if(!divise(e[9], p3)) {f = nudelta - 1; kod = 3; goto exit;}           
  1555.   /* III  */
  1556.       if (!divise(e[8], p3)) {f = nudelta - 2; kod = 4; goto exit;}          
  1557.   /* IV   */
  1558.       
  1559.       /* now for the last five cases... */
  1560.       
  1561.       if (!divise(e[5], p3))
  1562.         cumule(&v, &e, gun, gzero, gzero, p == 2 ? stoi(2) : modis(e[3], 9));
  1563.       /* p | a1, a2; p^2  | a3, a4; p^3 | a6 */
  1564.       a21 = aux((GEN)e[2], p, 1); a42 = aux((GEN)e[4], p, 2); 
  1565.       a63 = aux((GEN)e[5], p, 3);
  1566.       switch (numroots3(a21, a42, a63, p, &theroot))
  1567.         {
  1568.         case 3: f = nudelta - 4; kod = -1; goto exit;                        
  1569.  /* I0*  */
  1570.         case 2: /* calcul de nu */
  1571.           if (theroot) cumule(&v, &e, gun, stoi(theroot * p), gzero, gzero);
  1572.           /* p | a1; p^2  | a2, a3; p^3 | a4; p^4 | a6 */
  1573.           nu = 1;
  1574.           pk = p2;
  1575.           p2k = stoi(p * p * p * p);
  1576.           for(;;)
  1577.             {
  1578.               if (numroots2(1, aux2((GEN)e[3], p, pk),
  1579.                             -aux2((GEN)e[5], p, p2k), p, &theroot) == 2)
  1580.         break;
  1581.               if (theroot) cumule(&v, &e, gun, gzero, gzero, mulsi(theroot,pk));
  1582.               pk1 = pk; pk = mulsi(p, pk); p2k = mulsi(p, p2k);
  1583.               nu++;
  1584.               if (numroots2(a21, aux2((GEN)e[4], p, pk),
  1585.                             aux2((GEN)e[5], p, p2k), p, &theroot) == 2)
  1586.         break;
  1587.               if (theroot) cumule(&v, &e, gun, mulsi(theroot, pk1), gzero, gzero);
  1588.               p2k = mulsi(p, p2k);
  1589.               nu++;
  1590.             }
  1591.           f = nudelta - 4 - nu; kod = -4 - nu; goto exit;                    
  1592.   /* Inu* */
  1593.         case 1:
  1594.           if (theroot) cumule(&v, &e, gun, stoi(theroot * p), gzero, gzero); 
  1595.           /* p | a1; p^2  | a2, a3; p^3 | a4; p^4 | a6 */
  1596.           a32 = aux((GEN)e[3], p, 2); a64 = aux((GEN)e[5], p, 4);
  1597.           if(numroots2(1, a32, -a64, p, &theroot) == 2)
  1598.             {f = nudelta - 6; kod = -4; goto exit;}                          
  1599.   /* IV*  */
  1600.           if (theroot) cumule(&v, &e, gun, gzero, gzero, stoi(theroot*p*p));
  1601.           /* p | a1; p^2 | a2; p^3 | a3, a4; p^5 | a6 */
  1602.           p4 = mulii(p2, p2);
  1603.           if (!divise(e[4], p4)) {f = nudelta - 7; kod = -3; goto exit;}     
  1604.   /* III* */
  1605.           p6 = mulii(p4, p2);
  1606.           if (!divise(e[5], p6)) {f = nudelta - 8; kod = -2; goto exit;}     
  1607.   /* II*  */
  1608.           cumule(&v, &e, p1, gzero, gzero, gzero);  /* non minimal, on repart
  1609. pour un tour */
  1610.           nudelta -= 12;
  1611.         }
  1612.     }
  1613.  exit:
  1614.   tetpil = avma;
  1615.   result = cgetg(4, 17);
  1616.   result[1] = lstoi(f); result[2] = lstoi(kod); result[3] = lcopy(v);
  1617.   return gerepile(av, tetpil, result);
  1618. }
  1619.  
  1620. /*
  1621.   Calcul de toutes les fibres non elliptiques d'une courbe sur Z.
  1622.   Etant donne une courbe elliptique sous forme longue e, dont les coefficients
  1623.   sont entiers, renvoie une matrice dont les lignes sont de la forme
  1624.   [p, fp, kodp]. Il y a une ligne par diviseur premier du discriminant.
  1625.   Ceci n'est pas implemente dans Gp. Est-ce utile ?
  1626. */
  1627.  
  1628. GEN globaltatealgo(e)
  1629.     GEN  e;
  1630. {
  1631.   long k, l,av;
  1632.   GEN p1, p2, p3, prims, result;
  1633.  
  1634.   if ((typ(e) != 17) || (lg(e) < 14)) err(elliper1);
  1635.  
  1636.   prims = decomp(e[12]);
  1637.   l = lg(p1 = (GEN)prims[1]);
  1638.   p2 = (GEN)prims[2];
  1639.   if ((long)prims == avma) cgiv(prims);
  1640.   result = cgetg(4, 19);
  1641.   result[1] = (long)p1;
  1642.   result[2] = (long)p2;
  1643.   result[3] = (long)(p3 = cgetg(l, 18));
  1644.   for(k = 1; k < l; k++) p3[k] = lgeti(3);
  1645.   av = avma;
  1646.   for(k = 1; k < l; k++)
  1647.     {
  1648.       GEN q = localreduction(e, (GEN)p1[k]);
  1649.       affii(q[1], p2[k]);
  1650.       affii(q[2], p3[k]);
  1651.       avma = av;
  1652.     }
  1653.   return result;
  1654. }
  1655.  
  1656. /*
  1657.   Algorithme de reduction d'une courbe sur Q a sa forme standard.
  1658.   Etant donne une courbe elliptique sous forme longue e, dont les coefficients
  1659.   sont rationnels, renvoie son [N, [u, r, s, t]], ou N est le conducteur
  1660.   arithmetique de e, et [u, r, s, t] est le changement de variables qui reduit
  1661.   e a sa forme minimale globale dans laquelle a1 et a3 valent 0 ou 1, et a2
  1662.   vaut -1, 0 ou 1 et tel que u est un rationnel positif.
  1663. */
  1664.  
  1665. GEN globalreduction(e1)
  1666.     GEN e1;
  1667. {
  1668.   long i, k, l, m, tetpil, av = avma;
  1669.   GEN p1, prims, result, N = gun, u = gun, r, s, t;
  1670.   GEN v = cgetg(5, 17), a = cgetg(7, 17), e = cgetg(20, 17);
  1671.  
  1672.   if ((typ(e1) != 17) || (lg(e1) < 14)) err(elliper1);
  1673.  
  1674.   for(i = 1; i < 5; i++) a[i] = e1[i]; a[5] = zero; a[6] = e1[5];
  1675.   prims = decomp(denom(a));
  1676.   l = lg(p1 = (GEN)prims[1]);
  1677.   for(k = 1; k < l; k++)
  1678.     {
  1679.       int n = 0;
  1680.       for(i = 1; i < 7; i++)
  1681.         if (!gcmp0(a[i])) 
  1682.           {
  1683.             m = i * n + ggval(a[i], (GEN)p1[k]);
  1684.             while(m < 0) {n++; m += i;}
  1685.           }
  1686.       u = gmul(u, gpuigs(p1[k], n));
  1687.    }
  1688.   v[1] = linv(u); v[2] = v[3] = v[4] = zero;
  1689.   for(i = 1; i < 14; i++) e[i] = e1[i];
  1690.   for(; i < 20; i++) e[i] = zero;
  1691.   if (!gcmp1(u)) e = coordch(e, v);
  1692.   prims = decomp(e[12]);
  1693.   l = lg(p1 = (GEN)prims[1]);
  1694.   for(k = (signe(e[12]) < 0) + 1; k < l; k++)
  1695.     {
  1696.       GEN q = localreduction(e, (GEN)p1[k]);
  1697.       GEN v1 = (GEN)q[3];
  1698.       N = mulii(N, gpui(p1[k], q[1]));
  1699.       if (!gcmp1((GEN)v1[1])) cumule1(&v, &e, v1);
  1700.     }
  1701.   s = gdiventgs(e[1], -2);
  1702.   r = gdiventgs(gaddgs(gsub(gsub(e[2], gmul(s, e[1])), gsqr(s)), 1), -3);
  1703.   t = gdiventgs(gadd(e[3], gmul(r, e[1])), -2);
  1704.   cumule(&v, &e, gun, r, s, t);
  1705.   tetpil = avma;
  1706.   result = cgetg(3, 17); result[1] = lcopy(N); result[2] = lcopy(v);
  1707.   return gerepile(av, tetpil, result);
  1708. }
  1709.  
  1710. /* renvoie a_{k,l} avec les notations de Tate */
  1711.  
  1712. static int aux(ak, p, l)
  1713.     GEN ak;
  1714.     int p, l;
  1715. {
  1716.   long av = avma, pl = p, res;
  1717.   while(--l) pl *= p;
  1718.   res = itos(modis(divis(ak, pl), p));
  1719.   avma = av;
  1720.   return res;
  1721. }
  1722.  
  1723. static int aux2(ak, p, pl)
  1724.     GEN ak, pl;
  1725.     int p;
  1726. {
  1727.   long av = avma, res;
  1728.   res = itos(modis(divii(ak, pl), p));
  1729.   avma = av;
  1730.   return res;
  1731. }
  1732.  
  1733. /* renvoie le nombre de racines distinctes du polynome XXX + aXX + bX + c
  1734.    modulo p s'il y a une racine multiple, elle est renvoyee dans *mult */
  1735.  
  1736. static int numroots3(a, b, c, p, mult)
  1737.     int a, b, c, p, *mult;
  1738. {       
  1739.   if (p == 2)
  1740.     if ((c + a * b) & 1) return 3;
  1741.     else {*mult = b; return (a + b) & 1 ? 2 : 1;}
  1742.   else
  1743.     if (a % 3) {*mult = a * b; return (a * b * (1 - b) + c) % 3 ? 3 : 2;}
  1744.     else {*mult = -c; return b % 3 ? 3 : 1;}
  1745. }
  1746.  
  1747. /* idem pour aXX +bX + c */
  1748.  
  1749. static int numroots2(a, b, c, p, mult)
  1750.     int a, b, c, p, *mult;
  1751. {
  1752.   if (p == 2) {*mult = c; return b & 1 ? 2 : 1;}
  1753.   else {*mult = a * b; return (b * b - a * c) % 3 ? 2 : 1;}
  1754. }
  1755.  
  1756. /* cumule les effets de plusieurs chgts de variable. On traite a part les cas
  1757.    particuliers frequents, tels que soit u = 1, soit r' = s' = t' = 0 */
  1758.  
  1759. static void cumulev(vtotal, u, r, s, t)
  1760.     GEN *vtotal, u, r, s, t;
  1761. {
  1762.   long av = avma, tetpil;
  1763.   GEN temp, v = *vtotal, v3 = cgetg(5, 17);
  1764.   if (gcmp1(v[1]))
  1765.     {
  1766.       v3[1] = lcopy(u);
  1767.       v3[2] = ladd(v[2], r);
  1768.       v3[3] = ladd(v[3], s);
  1769.       av = avma;
  1770.       temp = gadd(v[4], gmul(v[3], r));
  1771.       tetpil = avma;
  1772.       v3[4] = lpile(av, tetpil, gadd(temp, t));
  1773.     }
  1774.   else if (gcmp0(r) && gcmp0(s) && gcmp0(t))
  1775.     {
  1776.       v3[1] = lmul(v[1], u);
  1777.       v3[2] = lcopy(v[2]);
  1778.       v3[3] = lcopy(v[3]);
  1779.       v3[4] = lcopy(v[4]);
  1780.     }
  1781.   else /* cas general */
  1782.     {
  1783.       v3[1] = lmul(v[1], u);
  1784.       temp = gmul(v[1],v[1]);
  1785.       v3[2] = ladd(gmul(temp, r), v[2]);
  1786.       v3[3] = ladd(gmul(v[1], s), v[3]);
  1787.       v3[4] = ladd(v[4], gmul(temp, gadd(gmul(v[1], t), gmul(v[3], r))));    
  1788.             
  1789.       tetpil = avma;
  1790.       v3 = gerepile(av, tetpil, gcopy(v3));
  1791.     }
  1792.   *vtotal = v3;
  1793. }
  1794.  
  1795. static void cumule(vtotal, e, u, r, s, t)
  1796.     GEN *vtotal, *e, u, r, s, t;
  1797. {
  1798.   long av = avma, tetpil;
  1799.   GEN v2 = cgetg(5, 17);
  1800.   v2[1] = (long)u; v2[2] = (long)r; v2[3] = (long)s; v2[4] = (long)t;
  1801.   tetpil = avma;
  1802.   *e = gerepile(av, tetpil, coordch(*e, v2));
  1803.   cumulev(vtotal, u, r, s, t);
  1804. }
  1805.  
  1806. static void cumule1(vtotal, e, v2)
  1807.     GEN *vtotal, *e, v2;
  1808. {
  1809.   *e = coordch(*e, v2);
  1810.   cumulev(vtotal, (GEN)v2[1], (GEN)v2[2], (GEN)v2[3], (GEN)v2[4]);
  1811. }
  1812.