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