home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Professional / OS2PRO194.ISO / os2 / progs / pari / pari_137 / src / polarit1.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-05-20  |  51.5 KB  |  1,915 lines
  1. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  2. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  3. /*~                                                                     ~*/
  4. /*~                       OPERATIONS ARITHMETIQUES                      ~*/
  5. /*~                                                                     ~*/
  6. /*~                           SUR LES POLYNOMES                         ~*/
  7. /*~                                                                     ~*/
  8. /*~                           (premiere partie)                         ~*/
  9. /*~                                                                     ~*/
  10. /*~                          copyright Babe Cool                        ~*/
  11. /*~                                                                     ~*/
  12. /*~                                                                     ~*/
  13. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  14. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  15.  
  16. # include "genpari.h"
  17.  
  18. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  19. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  20. /*                                                                 */
  21. /*                           DIVISIBILITE                          */
  22. /*                                                                 */
  23. /*                 Renvoie 1 si y divise x, 0 sinon .              */
  24. /*                                                                 */
  25. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  26. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  27.  
  28.  
  29. int gdivise(x,y)
  30.   
  31.    GEN   x,y;
  32.   
  33. {
  34.   long    avmacourant,i;
  35.   GEN   p1;
  36.  
  37.   avmacourant=avma;
  38.   p1=gmod(x,y);i=gcmp0(p1);
  39.   avma=avmacourant;
  40.   return i;
  41. }
  42.  
  43. int poldivis(x,y,z)
  44.      GEN x,y,*z;
  45. {
  46.   long av=avma;
  47.   GEN p1,p2;
  48.  
  49.   p1=poldivres(x,y,&p2);
  50.   if(signe(p2)) {avma=av;return 0;}
  51.   cgiv(p2);*z=p1;return 1;
  52. }
  53.  
  54.  
  55.  
  56. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  57. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  58. /*                                                                 */
  59. /*                            REDUCTION                            */
  60. /*                                                                 */
  61. /*        Met sous forme de fraction une nfraction; sous           */
  62. /*        forme de fr.rat une n.fr.rat; dans les autres cas        */
  63. /*        creation d'une copie .                                   */
  64. /*                                                                 */
  65. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  66. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  67.  
  68.  
  69.  
  70. GEN   gred(x)
  71.   
  72.    GEN   x;
  73.   
  74. {
  75.   long    tx,tetpil,l;
  76.   GEN   y,p1,x1,x2,x3;
  77.  
  78.   tx=typ(x);l=avma;
  79.   if ((tx==4)||(tx==5))
  80.   {
  81.     y=dvmdii(x[1],x[2],&p1);
  82.     if(!signe(p1)) cgiv(p1);
  83.     else
  84.     {
  85.       p1=mppgcd(x[2],p1);tetpil=avma;
  86.       y=cgetg(3,4);
  87.       if(gcmp1(p1))
  88.       {
  89.         y[1]=lcopy(x[1]);
  90.         y[2]=lcopy(x[2]);
  91.       }
  92.       else
  93.       {
  94.         y[1]=ldivii(x[1],p1);
  95.         y[2]=ldivii(x[2],p1);
  96.       }
  97.       y=gerepile(l,tetpil,y);
  98.     }
  99.   }
  100.   else if ((tx==13)||(tx==14))
  101.   {
  102.     x1=content(x[1]);x2=content(x[2]);x3=gdiv(x1,x2);
  103.     x1=(gcmp0(x1))?(GEN)x[1]:gdiv(x[1],x1);
  104.     x2=gdiv(x[2],x2);y=poldivres(x1,x2,&p1);
  105.     if(gcmp0(p1)) {tetpil=avma;return gerepile(l,tetpil,gmul(x3,y));}
  106.     else
  107.     {
  108.       p1=ggcd(x2,p1);y=cgetg(3,13);
  109.       if(isscalar(p1)) {y[1]=(long)x1;y[2]=(long)x2;}
  110.       else {y[1]=ldeuc(x1,p1);y[2]=ldeuc(x2,p1);}
  111.       x1=numer(x3);x2=denom(x3);tetpil=avma;
  112.       p1=cgetg(3,13);p1[1]=lmul(x1,y[1]);p1[2]=lmul(x2,y[2]);
  113.       return gerepile(l,tetpil,p1);
  114.     }
  115.   }
  116.   else y=gcopy(x);
  117.   return y;
  118. }
  119.  
  120. /*    REDUCTION SUR PLACE AVEC CHANGEMENT DE TYPE EVENTUEL  */
  121.  
  122. void gredsp(px)
  123.   
  124.    GEN   *px;
  125.   
  126. {
  127.   long    tx,l,l1;
  128.   GEN   x,y,p1;
  129.  
  130.   x= *px;tx=typ(x);
  131.   if ((tx==4)||(tx==5))
  132.   {
  133.     l=avma;
  134.     y=dvmdii(x[1],x[2],&p1);
  135.     if(!signe(p1))
  136.     {
  137.       cgiv(p1);l1=(long)(x+3);*px=gerepile(l1,l,y);
  138.     }
  139.     else
  140.     {
  141.       p1=mppgcd(x[2],p1);
  142.       if(!gcmp1(p1))
  143.       {
  144.         mpdivz(x[1],p1,x[1]);
  145.         mpdivz(x[2],p1,x[2]);
  146.       }
  147.       settyp(x,4);avma=l;
  148.     }
  149.   }
  150.   else if ((tx==13)||(tx==14))
  151.   {
  152.     l=avma;
  153.     y=poldivres(x[1],x[2],&p1);
  154.     if(!signe(p1))
  155.     {
  156.       cgiv(p1);l1=(long)(x+3);*px=gerepile(l1,l,y);
  157.     }
  158.     else
  159.     {
  160.       p1=primpart(ggcd(x[2],p1));
  161.       if(isnonscalar(p1))
  162.       {
  163.         gdeucz(x[1],p1,x[1]);
  164.         gdeucz(x[2],p1,x[2]);
  165.       }
  166.       settyp(y,13);avma=l;
  167.     }
  168.   }
  169. }
  170.  
  171.  
  172. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  173. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  174. /*                                                                 */
  175. /*               DIVISION EUCLIDIENNE DES POLYNOMES                */
  176. /*                                                                 */
  177. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  178. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  179.  
  180. GEN   gdeuc(x,y)
  181.   
  182.    GEN   x,y;
  183.   
  184. {
  185.   long  tx=typ(x),ty=typ(y),dx,dy,dz,vx,vy,i,j,l,tetpil,tetpil2,f;
  186.   GEN z,p1,p2,p3;
  187.  
  188.   if(ty<10) return gdiv(x,y);
  189.   if(tx<10) return gzero;
  190.   if((tx!=10)||(ty!=10)) err(poler1);
  191.   if(gcmp0(y)) err(poler4);
  192.   dx=lgef(x)-3;dy=lgef(y)-3;vx=varn(x);vy=varn(y);
  193.   if((vx>vy)||(dx<dy))
  194.   {z=cgetg(3,10);z[1]=2;z[2]=zero;setvarn(z,vy);}
  195.   else
  196.   {
  197.     if(vx<vy) z=gdiv(x,y);
  198.     else
  199.     {
  200.       dz=dx-dy;
  201.       z=cgetg(dz+3,10);z[1]=0x1000003+dz;setvarn(z,vx);
  202.       f=gcmp1(y[dy+2]);
  203.       if(f) z[dz+2]=lcopy(x[dx+2]);
  204.       else
  205.       z[dz+2]=ldiv(x[dx+2],y[dy+2]);
  206.       for(i=dx-1;i>=dy;--i)
  207.       {
  208.         l=avma;p1=(GEN)(x[i+2]);tetpil2=0;
  209.         for(j=i-dy+1;(j<=i)&&(j<=dz);j++)
  210.         {
  211.           p2=gmul(z[j+2],y[i-j+2]);
  212.           tetpil2=avma;p1=gsub(p1,p2);
  213.         }
  214.         tetpil=avma;
  215.         if(f)
  216.         {
  217.           if(tetpil2) z[i-dy+2]=lpile(l,tetpil2,p1);
  218.           else z[i-dy+2]=(long)p1;
  219.         }
  220.         else
  221.         {
  222.           p3=gdiv(p1,y[dy+2]);
  223.           z[i-dy+2]=lpile(l,tetpil,p3);
  224.         }
  225.       }
  226.     }
  227.   }
  228.   return z;
  229. }
  230.  
  231.  
  232. GEN   gres(x,y)
  233.   
  234.    GEN   x,y;
  235.   
  236. {
  237.   long  tx=typ(x),ty=typ(y),dx,dy,dz,i,j,k,f,l,l1,l2,tetpil,vx,vy;
  238.   GEN z,p1,p2,p3,p4;
  239.  
  240.   if(ty<10) return gzero;
  241.   if(tx<10) return gcopy(x);
  242.   if((tx!=10)||(ty!=10)) err(poler2);
  243.   if(gcmp0(y)) err(poler5);
  244.   dx=lgef(x)-3;dy=lgef(y)-3;vx=varn(x);vy=varn(y);
  245.   if((vx>vy)||(dx<dy)) z=gcopy(x);
  246.   else
  247.   {
  248.     if(vx<vy)
  249.     {z=cgetg(3,10);z[2]=zero;z[1]=2;setvarn(z,vx);}
  250.     else
  251.     {
  252.       dz=dx-dy;l1=avma;
  253.       p4=cgetg(dz+3,10);p4[1]=0x1000003+dz;setvarn(p4,vx);
  254.       p4[dz+2]=ldiv(x[dx+2],y[dy+2]);
  255.       for(i=dx-1;i>=dy;--i)
  256.       {
  257.         l=avma;p1=(GEN)(x[i+2]);
  258.         for(j=i-dy+1;(j<=i)&&(j<=dz);j++)
  259.         {
  260.           p2=gmul(p4[j+2],y[i-j+2]);
  261.           p1=gsub(p1,p2);
  262.         } tetpil=avma;p3=gdiv(p1,y[dy+2]);
  263.         p4[i-dy+2]=lpile(l,tetpil,p3);
  264.       }
  265.       l2=avma;f=1;
  266.       for(i=dy-1;(i>=0)&&f;--i)
  267.       {
  268.         l=avma;p1=(GEN)(x[i+2]);
  269.         for(j=0;(j<=i)&&(j<=dz);j++)
  270.         {
  271.           p2=gmul(p4[j+2],y[i-j+2]);
  272.           tetpil=avma;p1=gsub(p1,p2);
  273.         }
  274.         p1=gerepile(l,tetpil,p1);f=gcmp0(p1);
  275.       }
  276.       if(f)
  277.       {z=cgetg(3,10);z[1]=2;z[2]=zero;setvarn(z,vx);}
  278.       else
  279.       {
  280.         z=cgetg(i+4,10);z[1]=0x1000004+i;setvarn(z,vx);
  281.         z[i+3]=(long)p1;
  282.         for(k=i;k>=0;--k)
  283.         {
  284.           l=avma;p1=(GEN)(x[k+2]);
  285.           for(j=0;(j<=k)&&(j<=dz);j++)
  286.           {
  287.             p2=gmul(p4[j+2],y[k-j+2]);
  288.             tetpil=avma;p1=gsub(p1,p2);
  289.           }
  290.           z[k+2]=lpile(l,tetpil,p1);
  291.         }
  292.       }
  293.       z=gerepile(l1,l2,z);
  294.     }
  295.   }
  296.   return z;
  297. }
  298.  
  299.  
  300. GEN   poldivres(x,y,pr)
  301.   
  302.    GEN   x,y,*pr;
  303.   
  304. {
  305.   long  tx=typ(x),ty=typ(y),dx,dy,dz,i,j,k,f,l,tetpil,vx,vy;
  306.   GEN z,p1,p2,p3;
  307.  
  308.   if(ty<10) {*pr=gzero;return gdiv(x,y);}
  309.   if(tx<10) {*pr=gcopy(x);return gzero;}
  310.   if(tx!=10) err(poler3);
  311.   vx=varn(x);vy=gvar9(y);
  312.   if(vy>vx)
  313.   {
  314.     z=gdiv(x,y);p1=cgetg(3,10);*pr=p1;p1[1]=2;
  315.     p1[2]=zero;setvarn(p1,vx);
  316.   }
  317.   else
  318.   {
  319.     if(typ(y)!=10) err(poler3);
  320.     if(gcmp0(y)) err(poler6);
  321.     dx=lgef(x)-3;dy=lgef(y)-3;
  322.     if((vx>vy)||(dx<dy))
  323.     {
  324.       z=cgetg(3,10);z[1]=2;z[2]=zero;
  325.       setvarn(z,vy);*pr=gcopy(x);
  326.     }
  327.     else
  328.     {
  329.       dz=dx-dy;
  330.       z=cgetg(dz+3,10);z[1]=0x1000003+dz;setvarn(z,vx);
  331.       z[dz+2]=ldiv(x[dx+2],y[dy+2]);
  332.       for(i=dx-1;i>=dy;--i)
  333.       {
  334.         l=avma;p1=(GEN)(x[i+2]);
  335.         for(j=i-dy+1;(j<=i)&&(j<=dz);j++)
  336.         {
  337.           p2=gmul(z[j+2],y[i-j+2]);
  338.           p1=gsub(p1,p2);
  339.         }
  340.         tetpil=avma;p3=gdiv(p1,y[dy+2]);
  341.         z[i-dy+2]=lpile(l,tetpil,p3);
  342.       }
  343.       f=1;l=avma;
  344.       for(i=dy-1;(i>=0)&&f;--i)
  345.       {
  346.         l=avma;p1=(GEN)(x[i+2]);
  347.         for(j=0;(j<=i)&&(j<=dz);j++)
  348.         {
  349.           p2=gmul(z[j+2],y[i-j+2]);
  350.           tetpil=avma;p1=gsub(p1,p2);
  351.         } p1=gerepile(l,tetpil,p1);f=gcmp0(p1);
  352.       }
  353.       if(f)
  354.       {
  355.         avma=l;*pr=cgetg(3,10);(*pr)[1]=2;(*pr)[2]=zero;
  356.         setvarn(*pr,vx);
  357.       }
  358.       else
  359.       {
  360.         *pr=cgetg(i+4,10);(*pr)[1]=0x1000004+i;
  361.         setvarn(*pr,vx);(*pr)[i+3]=(long)p1;
  362.         for(k=i;k>=0;--k)
  363.         {
  364.           l=avma;p1=(GEN)(x[k+2]);
  365.           for(j=0;(j<=k)&&(j<=dz);j++)
  366.           {
  367.             p2=gmul(z[j+2],y[k-j+2]);
  368.             tetpil=avma;p1=gsub(p1,p2);
  369.           }
  370.           (*pr)[k+2]=lpile(l,tetpil,p1);
  371.         }
  372.       }
  373.     }
  374.   }
  375.   return z;
  376. }
  377.  
  378.  
  379.  
  380.  
  381. /*********************************************************************/
  382. /*********************************************************************/
  383. /*                                                                   */
  384. /*                         FONCTION BEZOUT                           */
  385. /*                                                                   */
  386. /*********************************************************************/
  387. /*********************************************************************/
  388.  
  389. GEN   bezoutpol(a,b,u,v)
  390.   
  391.    GEN    a,b,*u,*v;
  392.   
  393. {
  394.   GEN d,q,u1,u2,u3,w1,w2,w3,*bof;
  395.   long  av,av1,av2,va,dec;
  396.  
  397.   if ((typ(a)!=10)||(typ(b)!=10)) err(bezoutpoler);
  398.   if((va=varn(a))!=varn(b)) err(bezoutpoler);
  399.   if (lgef(b)>lgef(a))
  400.   {
  401.     u1=b;b=a;a=u1;bof=u;u=v;v=bof;
  402.   }
  403.   if (signe(b)) {
  404.   w1=a;w2=b;u1=polun[va];u2=gzero;
  405.   av=avma;
  406.   do
  407.     {
  408.       q=poldivres(w1,w2,&w3);
  409.       u3=gsub(u1,gmul(q,u2));
  410.       u1=u2;u2=u3;w1=w2;w2=w3;
  411.     }
  412.   while(signe(w3));
  413.   u3=gdiv(gsub(w1,gmul(u1,a)),b);av2=avma;
  414.   d=gdiv(w1,w3=leadingterm(w1));u2=gdiv(u1,w3);
  415.   u1=gdiv(u3,w3);av1=avma;dec=lpile(av,av2,0)>>2;
  416.   *u=adecaler(u2,av2,av1)?u2+dec:u2;*v=adecaler(u1,av2,av1)?u1+dec:u1;
  417.   if(adecaler(d,av2,av1)) d+=dec;
  418.   } /* fin du cas b non nul */
  419.   else {d=gcopy(a);*u= *v=polun[va];}
  420.   return d;
  421. }
  422.  
  423. GEN   polinvmod(x,y)
  424.   
  425.    GEN    x,y;
  426.   
  427. {
  428.   long  av,av1;
  429.   GEN u,v,d;
  430.  
  431.   av=avma;
  432.   d=bezoutpol(x,y,&u,&v);
  433.   if (gcmp0(d)||isnonscalar(d)) err(ginvmoder);
  434.   av1=avma;
  435.   return gerepile(av,av1,gdiv(u,d[2]));
  436. }
  437.  
  438.  
  439. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  440. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  441. /*                                                                 */
  442. /*                EVALUATION D'UN POLYNOME REEL                    */
  443. /*                                                                 */
  444. /*                                                                 */
  445. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  446. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  447.  
  448.  
  449. GEN   poleval(x,y)
  450.   
  451.    GEN   x,y;
  452.   
  453. {
  454.   long  l,av,tetpil,i,tx;
  455.   GEN z,p1,p2,p3,r,s;
  456.  
  457.   if((tx=typ(x))<10) return gcopy(x);
  458.   if(tx!=10) err(polevaler1);
  459.   l=lgef(x);
  460.   if (l==2) z=gzero;
  461.   else
  462.   {
  463.     if (l==3) z=gcopy(x[2]);
  464.     else
  465.     {
  466.       av=avma;p1=(GEN)x[l-1];
  467.       if(typ(y)!=6)
  468.       {
  469.         for (i=l-1;i>3;--i)
  470.         p1=gadd(gmul(p1,y),x[i-1]);
  471.         p1=gmul(y,p1);tetpil=avma;
  472.         z=gerepile(av,tetpil,gadd(p1,x[2]));
  473.       }
  474.       else
  475.       {
  476.         p2=(GEN)x[l-2];r=gadd(y[1],y[1]);
  477.         s=gnorm(y);
  478.         for(i=l-3;i>=2;--i)
  479.         {
  480.           p3=gadd(p2,gmul(r,p1));
  481.           p2=gsub(x[i],gmul(s,p1));
  482.           p1=p3;
  483.         }
  484.         p1=gmul(y,p1);tetpil=avma;
  485.         z=gerepile(av,tetpil,gadd(p1,p2));
  486.       }
  487.     }
  488.   }
  489.   return z;
  490. }
  491.  
  492.  
  493.  
  494. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  495. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  496. /*                                                                 */
  497. /*                         RACINES COMPLEXES                       */
  498. /*                                                                 */
  499. /*        l represente la longueur voulue pour les parties         */
  500. /*            reelles et imaginaires des racines de x              */
  501. /*                                                                 */
  502. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  503. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  504.  
  505. GEN   roots(x,l)
  506.   
  507.    long  l;
  508.    GEN   x;
  509.   
  510. {
  511.   long  av,i,j,f,g,fr,deg,l0,l1,l2,l3,l4,ln,dec;
  512.   long  exc,expmin,m,deg0,k,ti,h,ii,e,e1,emax,v;
  513.   GEN y,xc,xd0,xd,xdabs,p1,p2,p3,p4,p5,p6,p7,p8;
  514.   GEN p9,p10,p11,p12,p14,p15,pa,pax,pb,pp,pq,ps;
  515.  
  516.   if(typ(x)!=10) err(poler7);
  517.   else
  518.   {
  519.     v=varn(x);deg0=lgef(x)-3;expmin=(2-l)*32+12;
  520.     if(!signe(x)) err(poler8);
  521.     y=cgetg(deg0+1,18);if(!deg0) return y;
  522.     for(i=1;i<=deg0;i++)
  523.     {
  524.       p1=cgetg(3,6);p1[1]=lgetr(l);
  525.       p1[2]=lgetr(l);y[i]=(long)p1;
  526.     }
  527.     g=1;f=1;
  528.     for(i=2;i<=deg0+2;i++)
  529.     {
  530.       ti=typ(x[i]);
  531.       if(ti==8)
  532.       {
  533.         p2=(GEN)((GEN)((GEN)x[i])[1])[2];
  534.         if(gsigne(p2)>0) g=0;
  535.       }
  536.       else
  537.       if(ti>5) g=0;
  538.     }
  539.     l1=avma;p2=cgetg(3,6);
  540.     p2[1]=lmppi(3);
  541.     p2[2]=ldivrs(p2[1],10);
  542.     p11=cgetg(4,10);p11[1]=0x1000004;setvarn(p11,v);p11[3]=un;
  543.     p12=cgetg(5,10);p12[1]=0x1000005;setvarn(p12,v);p12[4]=un;
  544.     for(i=2;(i<=deg0+2)&&(gcmp0(x[i]));i++)
  545.     gaffsg(0,y[i-1]);k=i-2;
  546.     if(k!=deg0)
  547.     {
  548.       if(k)
  549.       {
  550.         j=deg0+3-k;pax=cgetg(j,10);pax[1]=0x1000000+j;
  551.         setvarn(pax,v);
  552.         for(i=2;i<j;i++)
  553.         pax[i]=x[i+k];
  554.       }
  555.       else pax=x;
  556.       xd0=deriv(pax,v);pp=ggcd(pax,xd0);m=1;pa=pax;
  557.       h=isnonscalar(pp);if(h) pq=gdeuc(pax,pp);
  558.       do
  559.       {
  560.         if(h)
  561.         {
  562.           pa=pp;pb=pq;
  563.           pp=ggcd(pa,deriv(pa,v));h=isnonscalar(pp);
  564.           if(h) pq=gdeuc(pa,pp);else pq=pa;
  565.           ps=gdeuc(pb,pq);
  566.         }
  567.         else ps=pa;
  568.         /* calcul des racines d'ordre exactement m */
  569.         deg=lgef(ps)-3;
  570.         if(deg)
  571.         {
  572.           l3=avma;e=gexpo(ps[deg+2]);emax=e;
  573.           for(i=2;i<deg+2;i++)
  574.           {
  575.             p3=(GEN)(ps[i]);
  576.             if(!gcmp0(p3))
  577.             {
  578.               e1=gexpo(p3);
  579.               if(e1>emax) emax=e1;
  580.             }
  581.           }
  582.           e=emax-e;if(e<0) e=0;avma=l3;
  583.           if(ps!=pax) xd0=deriv(ps,v);
  584.           xdabs=cgetg(deg+2,10);xdabs[1]=xd0[1];
  585.           for(i=2;i<deg+2;i++)
  586.           {
  587.             l3=avma;p3=(GEN)xd0[i];p4=gabs(greal(p3));
  588.             p5=gabs(gimag(p3),l);l4=avma;
  589.             xdabs[i]=lpile(l3,l4,gadd(p4,p5));
  590.           }
  591.           l0=avma;xc=gcopy(ps);xd=gcopy(xd0);l2=avma;
  592.           for(i=1;i<=deg;i++)
  593.           {
  594.             if(i==deg)
  595.             {
  596.               p1=(GEN)y[k+m*i];
  597.               gdivz(gneg(xc[2]),xc[3],p1);
  598.               p14=(GEN)(p1[1]);p15=(GEN)(p1[2]);
  599.             }
  600.             else
  601.             {
  602.               p3=gshift(p2,e);p4=poleval(xc,p3);
  603.               p5=gnorm(p4);exc=0;
  604.               while(exc>= -20)
  605.               {
  606.                 p6=poleval(xd,p3);p7=gneg(gdiv(p4,p6));
  607.                 f=1;l3=avma;if(gcmp0(p5)) exc= -32;
  608.                 else exc=expo(gnorm(p7))-expo(gnorm(p3));
  609.                 avma=l3;
  610.                 for(j=1;(j<=10)&&f;j++)
  611.                 {
  612.                   p8=gadd(p3,p7);p9=poleval(xc,p8);
  613.                   p10=gnorm(p9);
  614.                   f=(cmprr(p10,p5)>=0)&&(exc>= -20);
  615.                   if(f) {gshiftz(p7,-2,p7);avma=l3;}
  616.                 }
  617.                 if(f) err(poler9);
  618.                 else
  619.                 {
  620.                   av=avma;dec=lpile(l2,l3,0)>>2;
  621.                   p3=adecaler(p8,l3,av)?p8+dec:p8;
  622.           p4=adecaler(p9,l3,av)?p9+dec:p9;
  623.           p5=adecaler(p10,l3,av)?p10+dec:p10;
  624.                 }
  625.               }
  626.               p1=(GEN)y[k+m*i];setlg(p1[1],3);
  627.               setlg(p1[2],3);gaffect(p3,p1);avma=l2;
  628.               p14=(GEN)(p1[1]);p15=(GEN)(p1[2]);
  629.               for(ln=4;ln<=l;ln=(ln<<1)-2)
  630.               {
  631.                 setlg(p14,ln);setlg(p15,ln);
  632.                 if(gcmp0(p14))
  633.                 {settyp(p14,1);p14[1]=2;}
  634.                 if(gcmp0(p15))
  635.                 {settyp(p15,1);p15[1]=2;}
  636.                 p4=poleval(xc,p1);p5=poleval(xd,p1);
  637.                 p6=gneg(gdiv(p4,p5));
  638.                 settyp(p14,2);settyp(p15,2);
  639.                 gaffect(gadd(p1,p6),p1);avma=l2;
  640.               }
  641.             }
  642.             setlg(p14,l);setlg(p15,l);
  643.             p7=gcopy(p1);
  644.             p14=(GEN)(p7[1]);p15=(GEN)(p7[2]);
  645.             setlg(p14,l+1);setlg(p15,l+1);
  646.             if(gcmp0(p14))
  647.             {settyp(p14,1);p14[1]=2;}
  648.             if(gcmp0(p15))
  649.             {settyp(p15,1);p15[1]=2;}
  650.             for(ii=1;ii<=5;ii++)
  651.             {
  652.               p4=poleval(ps,p7);p5=poleval(xd0,p7);
  653.               p6=gneg(gdiv(p4,p5));p7=gadd(p7,p6);
  654.               p14=(GEN)(p7[1]);p15=(GEN)(p7[2]);
  655.               if(gcmp0(p14))
  656.               {settyp(p14,1);p14[1]=2;}
  657.               if(gcmp0(p15))
  658.               {settyp(p15,1);p15[1]=2;}
  659.             }
  660.             gaffect(p7,p1);p6=gdiv(p4,poleval(xdabs,gabs(p7,l)));
  661.             avma=l2;
  662.             /*          if(!gcmp0(p6))
  663.                       if(gexpo(p6)>=expmin) err(poler10);*/
  664.             if((expo(p1[2])<expmin)&&g)
  665.             {
  666.               gaffect(zero,p1[2]);
  667.               for(j=1;j<m;j++)
  668.               gaffect(p1,y[k+(i-1)*m+j]);
  669.               p11[2]=lneg(p1[1]);l4=avma;
  670.               xc=gerepile(l0,l4,gdeuc(xc,p11));
  671.             }
  672.             else
  673.             {
  674.               for(j=1;j<m;j++)
  675.               gaffect(p1,y[k+(i-1)*m+j]);
  676.               if(g)
  677.               {
  678.                 p1=gconj(p1);
  679.                 for(j=1;j<=m;j++)
  680.                 gaffect(p1,y[k+i*m+j]);i++;
  681.                 p12[2]=lnorm(p1);
  682.                 p12[3]=lmulsg(-2,p1[1]);
  683.                 l4=avma;
  684.                 xc=gerepile(l0,l4,gdeuc(xc,p12));
  685.               }
  686.               else
  687.               {
  688.                 p11[2]=lneg(p1);l4=avma;
  689.                 xc=gerepile(l0,l4,gdeuc(xc,p11));
  690.               }
  691.             }
  692.             xd=deriv(xc,v);l2=avma;
  693.           }
  694.           k=k+deg*m;
  695.         }
  696.         m++;
  697.       }
  698.       while (k!=deg0);
  699.     }
  700.     avma=l1;
  701.     if(g&&(deg0>1))
  702.     {
  703.       for(j=2;j<=deg0;j++)
  704.       {
  705.         p1=(GEN)y[j];fr=gcmp0(p1[2]);
  706.         i=j-1;
  707.         if(fr)
  708.         {
  709.           p2=(GEN)y[i];f=!gcmp0(p2[2]);
  710.           if(!f) f=(cmprr(p2[1],p1[1])>0);
  711.           for(;(i>0)&&f;--i)
  712.           {
  713.             y[i+1]=y[i];
  714.             p2=(GEN)y[i];f=!gcmp0(p2[2]);
  715.             if(!f) f=(cmprr(p2[1],p1[1])>0);
  716.           }
  717.         }
  718.         y[i+1]=(long)p1;
  719.       }
  720.     }
  721.   }
  722.   return y;
  723. }
  724.  
  725. GEN   rootslong(x,l)
  726.   
  727.    long  l;
  728.    GEN   x;
  729.   
  730. {
  731.   long  av,i,j,f,g,fr,deg,l0,l1,l2,l3,l4,ln,dec;
  732.   long  exc,expmin,m,deg0,k,ti,h,ii,e,e1,emax,v;
  733.   GEN y,xc,xd0,xd,xdabs,p1,p2,p3,p4,p5,p6,p7,p8;
  734.   GEN p9,p10,p11,p12,p14,p15,pa,pax,pb,pp,pq,ps;
  735.  
  736.   if(typ(x)!=10) err(poler7);
  737.   else
  738.   {
  739.     v=varn(x);deg0=lgef(x)-3;expmin=(2-l)*32+12;
  740.     if(!signe(x)) err(poler8);
  741.     y=cgetg(deg0+1,18);if(!deg0) return y;
  742.     for(i=1;i<=deg0;i++)
  743.     {
  744.       p1=cgetg(3,6);p1[1]=lgetr(l);
  745.       p1[2]=lgetr(l);y[i]=(long)p1;
  746.     }
  747.     g=1;f=1;
  748.     for(i=2;i<=deg0+2;i++)
  749.     {
  750.       ti=typ(x[i]);
  751.       if(ti==8)
  752.       {
  753.         p2=(GEN)((GEN)((GEN)x[i])[1])[2];
  754.         if(gcmpgs(p2,0)>0) g=0;
  755.       }
  756.       else
  757.       if(ti>5) g=0;
  758.     }
  759.     l1=avma;p2=cgetg(3,6);
  760.     p2[1]=lmppi(l);
  761.     p2[2]=ldivrs(p2[1],10);
  762.     p11=cgetg(4,10);p11[1]=0x1000004;setvarn(p11,v);p11[3]=un;
  763.     p12=cgetg(5,10);p12[1]=0x1000005;setvarn(p12,v);p12[4]=un;
  764.     for(i=2;(i<=deg0+2)&&(gcmp0(x[i]));i++)
  765.     gaffsg(0,y[i-1]);k=i-2;
  766.     if(k!=deg0)
  767.     {
  768.       if(k)
  769.       {
  770.         j=deg0+3-k;pax=cgetg(j,10);pax[1]=0x1000000+j;
  771.         setvarn(pax,v);
  772.         for(i=2;i<j;i++)
  773.         pax[i]=x[i+k];
  774.       }
  775.       else pax=x;
  776.       xd0=deriv(pax,v);pp=ggcd(pax,xd0);m=1;pa=pax;
  777.       h=isnonscalar(pp);if(h) pq=gdeuc(pax,pp);
  778.       do
  779.       {
  780.         if(h)
  781.         {
  782.           pa=pp;pb=pq;
  783.           pp=ggcd(pa,deriv(pa,v));h=isnonscalar(pp);
  784.           if(h) pq=gdeuc(pa,pp);else pq=pa;
  785.           ps=gdeuc(pb,pq);
  786.         }
  787.         else ps=pa;
  788.         /* calcul des racines d'ordre exactement m */
  789.         deg=lgef(ps)-3;
  790.         if(deg)
  791.         {
  792.           l3=avma;e=gexpo(ps[deg+2]);emax=e;
  793.           for(i=2;i<deg+2;i++)
  794.           {
  795.             p3=(GEN)(ps[i]);
  796.             if(!gcmp0(p3))
  797.             {
  798.               e1=gexpo(p3);
  799.               if(e1>emax) emax=e1;
  800.             }
  801.           }
  802.           e=emax-e;if(e<0) e=0;avma=l3;
  803.           if(ps!=pax) xd0=deriv(ps,v);
  804.           xdabs=cgetg(deg+2,10);xdabs[1]=xd0[1];
  805.           for(i=2;i<deg+2;i++)
  806.           {
  807.             l3=avma;p3=(GEN)xd0[i];p4=gabs(greal(p3));
  808.             p5=gabs(gimag(p3),l);l4=avma;
  809.             xdabs[i]=lpile(l3,l4,gadd(p4,p5));
  810.           }
  811.           l0=avma;xc=gcopy(ps);xd=gcopy(xd0);l2=avma;
  812.           for(i=1;i<=deg;i++)
  813.           {
  814.             if(i==deg)
  815.             {
  816.               p1=(GEN)y[k+m*i];
  817.               gdivz(gneg(xc[2]),xc[3],p1);
  818.               p14=(GEN)(p1[1]);p15=(GEN)(p1[2]);
  819.             }
  820.             else
  821.             {
  822.               p3=gshift(p2,e);p4=poleval(xc,p3);
  823.               p5=gnorm(p4);exc=0;
  824.               while(exc>= -20)
  825.               {
  826.                 p6=poleval(xd,p3);p7=gneg(gdiv(p4,p6));
  827.                 f=1;l3=avma;if(gcmp0(p5)) exc= -32;
  828.                 else exc=expo(gnorm(p7))-expo(gnorm(p3));
  829.                 avma=l3;
  830.                 for(j=1;(j<=50)&&f;j++)
  831.                 {
  832.                   p8=gadd(p3,p7);p9=poleval(xc,p8);
  833.                   p10=gnorm(p9);
  834.                   f=(cmprr(p10,p5)>=0)&&(exc>= -20);
  835.                   if(f) {gshiftz(p7,-2,p7);avma=l3;}
  836.                 }
  837.                 if(f) err(poler9);
  838.                 else
  839.                 {
  840.                   av=avma;dec=lpile(l2,l3,0)>>2;
  841.                   p3=adecaler(p8,l3,av)?p8+dec:p8;
  842.           p4=adecaler(p9,l3,av)?p9+dec:p9;
  843.           p5=adecaler(p10,l3,av)?p10+dec:p10;
  844.                 }
  845.               }
  846.               p1=(GEN)y[k+m*i];gaffect(p3,p1);avma=l2;
  847.               p14=(GEN)(p1[1]);p15=(GEN)(p1[2]);
  848.               for(ln=4;ln<=l;ln=(ln<<1)-2)
  849.               {
  850.                 if(gcmp0(p14))
  851.                 {settyp(p14,1);p14[1]=2;}
  852.                 if(gcmp0(p15))
  853.                 {settyp(p15,1);p15[1]=2;}
  854.                 p4=poleval(xc,p1);p5=poleval(xd,p1);
  855.                 p6=gneg(gdiv(p4,p5));
  856.                 settyp(p14,2);settyp(p15,2);
  857.                 gaffect(gadd(p1,p6),p1);avma=l2;
  858.               }
  859.             }
  860.             p7=gcopy(p1);
  861.             p14=(GEN)(p7[1]);p15=(GEN)(p7[2]);
  862.             setlg(p14,l+1);setlg(p15,l+1);
  863.             if(gcmp0(p14))
  864.             {settyp(p14,1);p14[1]=2;}
  865.             if(gcmp0(p15))
  866.             {settyp(p15,1);p15[1]=2;}
  867.             for(ii=1;ii<=((e<<5)+2);ii<<=1)
  868.             {
  869.               p4=poleval(ps,p7);p5=poleval(xd0,p7);
  870.               p6=gneg(gdiv(p4,p5));p7=gadd(p7,p6);
  871.               p14=(GEN)(p7[1]);p15=(GEN)(p7[2]);
  872.               if(gcmp0(p14))
  873.               {settyp(p14,1);p14[1]=2;}
  874.               if(gcmp0(p15))
  875.               {settyp(p15,1);p15[1]=2;}
  876.             }
  877.             gaffect(p7,p1);p6=gdiv(p4,poleval(xdabs,gabs(p7,l)));
  878.             avma=l2;
  879.             /*          if(!gcmp0(p6))
  880.                       if(gexpo(p6)>=expmin) err(poler10);*/
  881.             if((expo(p1[2])<expmin)&&g)
  882.             {
  883.               gaffect(zero,p1[2]);
  884.               for(j=1;j<m;j++)
  885.               gaffect(p1,y[k+(i-1)*m+j]);
  886.               p11[2]=lneg(p1[1]);l4=avma;
  887.               xc=gerepile(l0,l4,gdeuc(xc,p11));
  888.             }
  889.             else
  890.             {
  891.               for(j=1;j<m;j++)
  892.               gaffect(p1,y[k+(i-1)*m+j]);
  893.               if(g)
  894.               {
  895.                 p1=gconj(p1);
  896.                 for(j=1;j<=m;j++)
  897.                 gaffect(p1,y[k+i*m+j]);i++;
  898.                 p12[2]=lnorm(p1);
  899.                 p12[3]=lmulsg(-2,p1[1]);
  900.                 l4=avma;
  901.                 xc=gerepile(l0,l4,gdeuc(xc,p12));
  902.               }
  903.               else
  904.               {
  905.                 p11[2]=lneg(p1);l4=avma;
  906.                 xc=gerepile(l0,l4,gdeuc(xc,p11));
  907.               }
  908.             }
  909.             xd=deriv(xc,v);l2=avma;
  910.           }
  911.           k=k+deg*m;
  912.         }
  913.         m++;
  914.       }
  915.       while (k!=deg0);
  916.     }
  917.     avma=l1;
  918.     if(g&&(deg0>1))
  919.     {
  920.       for(j=2;j<=deg0;j++)
  921.       {
  922.         p1=(GEN)y[j];fr=gcmp0(p1[2]);
  923.         i=j-1;
  924.         if(fr)
  925.         {
  926.           p2=(GEN)y[i];f=!gcmp0(p2[2]);
  927.           if(!f) f=(cmprr(p2[1],p1[1])>0);
  928.           for(;(i>0)&&f;--i)
  929.           {
  930.             y[i+1]=y[i];
  931.             p2=(GEN)y[i];f=!gcmp0(p2[2]);
  932.             if(!f) f=(cmprr(p2[1],p1[1])>0);
  933.           }
  934.         }
  935.         y[i+1]=(long)p1;
  936.       }
  937.     }
  938.   }
  939.   return y;
  940. }
  941.  
  942.  
  943. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  944. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  945. /*                                                                 */
  946. /*      Recherche des racines modulo p ( par verif   f(x)=0 )      */
  947. /*                                                                 */
  948. /*      (retourne le vecteur horizontal dont les composantes       */
  949. /*       sont les racines (eventuellement vecteur a 0 comp.)       */
  950. /*                                                                 */
  951. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  952. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  953.  
  954. GEN     rootmod2(f,p)
  955.  
  956.   GEN     f,p;
  957.  
  958. {
  959.   GEN g,y,z,ss;
  960.   long vf,av,av1,av2,deg,s,nbrac,pasfini;
  961.  
  962.   if((typ(f)!=10)||gcmp0(f)) err(factmoder);
  963.   y=(GEN)(f[2]);vf=varn(f);
  964.   if((typ(y)!=3)||!gegal(y[1],p))
  965.     {p=gcopy(p);av=avma;f=gmul(f,gmodulcp(un,p));}
  966.   else 
  967.     av=avma;
  968.   deg=lgef(f)-3;
  969.   s=0;nbrac=0;
  970.   y=cgetg(deg+1,17);
  971.   av1=avma;
  972.   do
  973.     {
  974.       if(av1==avma)
  975.     {
  976.       ss=stoi(s);
  977.       pasfini=(gcmp(p,ss)>0);
  978.     }
  979.       else av1=avma;
  980.       av2=avma;
  981.       z=poleval(f,ss);
  982.       if(gcmp0(z))
  983.     {
  984.       avma=av2;
  985.       nbrac++;
  986.       y[nbrac]=(long)gmodulcp(ss,p);
  987.       f=gdiv(f,gsub(polx[vf],ss));
  988.     }
  989.       else
  990.     {
  991.       avma=av1;s++;
  992.     }
  993.     }
  994.   while((nbrac<deg-1)&&pasfini);
  995.   if (!nbrac) {avma=av;return cgetg(1,17);}
  996.   if (nbrac==(deg-1))
  997.     {
  998.       nbrac++;y[nbrac]=lneg(gdiv(f[2],f[3]));
  999.     }
  1000.   g=cgetg(nbrac+1,17);
  1001.   for(s=1;s<=nbrac;s++) g[s]=y[s];
  1002.   av1=avma;return gerepile(av,av1,gcopy(g));
  1003. }
  1004.  
  1005. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1006. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1007. /*                                                                 */
  1008. /*          Recherche intelligente des racines modulo p            */
  1009. /*                                                                 */
  1010. /*      (retourne le vecteur horizontal dont les composantes       */
  1011. /*       sont les racines (eventuellement vecteur a 0 comp.)       */
  1012. /*                                                                 */
  1013. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1014. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1015.  
  1016.  
  1017. GEN rootmod(f,p)
  1018.      GEN f,p;
  1019.  
  1020. {
  1021.   GEN y,unmodp,pol,xun,g,a,b,d,e,u,h,q,p1;
  1022.   long av,tetpil,vf,n,i,j,la,flag,lb,rac[10];
  1023.  
  1024.   if((typ(f)!=10)||gcmp0(f)) err(factmoder);
  1025.   y=(GEN)(f[2]);vf=varn(f);
  1026.   if((typ(y)!=3)||!gegal(y[1],p))
  1027.     {p=gcopy(p);av=avma;f=gmul(f,gmodulcp(un,p));}
  1028.   else av=avma;
  1029.   unmodp=gmodulcp(gun,p);
  1030.   if(gegal(p,deux))
  1031.     {
  1032.       j=0;if(gcmp0(f[2])) j++;
  1033.       if(gcmp0(gsubst(f,vf,unmodp))) j+=2;
  1034.       avma=av;
  1035.       switch(j)
  1036.     {
  1037.     case 0: y=cgetg(1,17);break;
  1038.     case 1: y=cgetg(2,17);y[1]=lmodulcp(gzero,p);break;
  1039.     case 2: y=cgetg(2,17);y[1]=lmodulcp(gun,p);break;
  1040.     case 3: y=cgetg(3,17);y[1]=lmodulcp(gzero,p);
  1041.       y[2]=lmodulcp(gun,p);
  1042.     }
  1043.       return y;
  1044.     }
  1045.   if(!gcmpgs(p,4))
  1046.     {
  1047.       j=0;if(gcmp0(f[2])) {j++;rac[j]=0;}
  1048.       p1=unmodp;
  1049.       for(i=1;i<=3;i++)
  1050.     {
  1051.       if(gcmp0(gsubst(f,vf,p1))) {j++;rac[j]=i;}
  1052.       if(i<3) p1=gadd(p1,unmodp);
  1053.     }
  1054.       avma=av;y=cgetg(j+1,17);
  1055.       for(i=1;i<=j;i++) y[i]=lmodulcp(stoi(rac[i]),p);
  1056.       return y;
  1057.     }
  1058.   pol=gmul(unmodp,polx[vf]);
  1059.   xun=gmodulcp(pol,f);g=ggcd((gsub(gpui(xun,p),xun))[2],f);
  1060.   n=lgef(g)-3;if(!n) {avma=av;y=cgetg(1,17);return y;}
  1061.   y=cgetg(n+1,17);
  1062.   if(gcmp0(g[2]))
  1063.     {
  1064.       y[1]=zero;g=gdiv(g,polx[vf]);
  1065.       if(lgef(g)>3) y[2]=(long)g;
  1066.       j=2;
  1067.     }
  1068.   else {y[1]=(long)g;j=1;}
  1069.   while(j<=n)
  1070.     {
  1071.       a=(GEN)y[j];la=lgef(a)-3;
  1072.       if(la==1)
  1073.     {
  1074.       y[j]=(gneg(gdiv(a[2],a[3])))[2];j++;
  1075.     }
  1076.       else if(la==2)
  1077.     {
  1078.       d=gsub(gmul(a[3],a[3]),gmul2n(gmul(a[2],a[4]),2));
  1079.       e=gsqrt(d);u=gdiv(gun,gmul2n(a[4],1));
  1080.       y[j]=(gmul(u,gsub(e,a[3])))[2];y[j+1]=(gmul(u,gneg(gadd(e,a[3]))))[2];
  1081.       j+=2;
  1082.     }
  1083.       else
  1084.     {
  1085.       flag=1;h=pol;q=shifti(subis(p,1),-1);
  1086.       while(flag)
  1087.         {
  1088.           p1=gmodulcp(h,a);b=ggcd((gsub(gpui(p1,q),gun))[2],a);
  1089.           lb=lgef(b)-3;
  1090.           if(lb&&(lb<la))
  1091.         {
  1092.           flag=0;y[j]=(long)b;y[j+lb]=ldiv(a,b);
  1093.         }
  1094.           else h=gadd(h,unmodp);
  1095.         }
  1096.     }
  1097.     }
  1098.   y=sort(y);tetpil=avma;return gerepile(av,tetpil,gmul(unmodp,y));
  1099. }
  1100.  
  1101. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1102. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1103. /*                                                                 */
  1104. /*                     FACTORISATION MODULO p                      */
  1105. /*                                                                 */
  1106. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1107. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1108.  
  1109. GEN     factmod(f,p)
  1110.      
  1111.      GEN     f,p;
  1112.      
  1113. {
  1114.   long  i,j,k,d,e,vf,expos[100],psim,nbfact,nbf,av,pk,tetpil,calc;
  1115.   GEN   y,t[100],f1,f2,f3,df1,df2,g,g1,g2;
  1116.   GEN   xmod,u,v,pd,q,a0;
  1117.   
  1118.   if((typ(f)!=10)||gcmp0(f)) err(factmoder);
  1119.   if((lgef(p)>3)||((lgef(p)==3)&&(cmpis(p,VERYBIGINT)>0))) 
  1120.     err(impl,"factmod for primes >2^31");
  1121.   a0=(GEN)(f[2]);vf=varn(f);
  1122.   if((typ(a0)!=3)||!gegal(a0[1],p))
  1123.     {p=gcopy(p);av=avma;f=gmul(f,gmodulcp(un,p));}
  1124.   else av=avma;
  1125.   if(lgef(f)==3)
  1126.     {avma=av;y=cgetg(3,19);y[1]=lgetg(1,18);y[2]=lgetg(1,18);return y;}
  1127.   
  1128.   /*  1ere etape : trouver les facteurs square-free produit
  1129.       des facteurs premiers de meme multiplicite    */
  1130.  
  1131.   f1=f;e=0;nbfact=1;pk=1;df1=deriv(f1,vf);calc=0;
  1132.  
  1133.   do
  1134.     {
  1135. /*printf("debut 1re etape\n");*/
  1136.       e+=pk;
  1137.       while (gcmp0(df1))
  1138.         {
  1139.           pk *= (psim=itos(p));e=pk;
  1140.           j=(lgef(f1)-3)/psim+3;f2=cgetg(j,10);
  1141.           f2[1]=0x1000000+j;setvarn(f2,vf);
  1142.           for(i=2;i<j;i++)
  1143.             f2[i]=lcopy(f1[psim*(i-2)+2]);
  1144.           f1=f2;df1=deriv(f1,vf);calc=0;
  1145.     }
  1146.       if(calc) f2=f3;else f2=ggcd(f1,df1);calc=1;
  1147.       if (lgef(f2)<4) u=f1;
  1148.       else
  1149.         {
  1150.           g1=gdiv(f1,f2);
  1151.           if (gcmp0(df2=deriv(f2,vf)))
  1152.             {
  1153.               u=g1;f3=f2;
  1154.         }
  1155.           else
  1156.             {
  1157.               f3=ggcd(f2,df2);
  1158.               if (lgef(f3)<4) u=gdiv(g1,f2);
  1159.               else
  1160.                 {
  1161.                   g2=gdiv(f2,f3);u=gdiv(g1,g2);
  1162.          }
  1163.         }
  1164.     }
  1165.       
  1166.       /*  2eme etape : 
  1167.           Ici u est un polynome square-free
  1168.           (produit des facteurs premiers de meme multiplicite  e :
  1169.           trouver les facteurs (square-free) produit des facteurs de meme
  1170.       degre d */
  1171.       
  1172.       d=0;pd=gun;
  1173.       xmod=gmodulcp(polx[vf],u);v=xmod;
  1174.       while(d<(lgef(u)-3)>>1)
  1175.         {
  1176.                                        /*printf("debut 2me etape\n");*/
  1177.           d++;
  1178.           pd=mulii(pd,p);
  1179.           q=shifti(subis(pd,1),-1);
  1180.           v=gpui(v,p);
  1181.           g=ggcd((gsub(v,xmod))[2],u);
  1182.       
  1183.           if (lgef(g)>3)
  1184.             {
  1185.                                       /*printf("debut 3me etape\n");*/
  1186.           
  1187. /*  3eme etape :
  1188. Ici g est produit de pol irreductibles ayant tous le meme degre d;*/
  1189.  
  1190.           t[nbfact]=g;j=nbfact+(lgef(g)-3)/d;
  1191.           split(p[2],t+nbfact,d,p[2],q);
  1192. /* le premier parametre est un entier variable m qui sera converti en un
  1193.    polynome w dont les coeff sont ses digits en base p (initialement m = p --> X)
  1194.    pour faire pgcd de g avec w^(p^d-1)/2 jusqu'a casser. */
  1195.           for(;nbfact<j;expos[nbfact++]=e);
  1196.           u=gdiv(u,g);v=gmodulcp(v[2],u);
  1197.         }
  1198.                                           /*printf("fin 3me etape\n");*/
  1199.  
  1200.     }                                 /*printf("fin 2me etape\n");*/
  1201.       if (lgef(u)>3) {t[nbfact]=u;expos[nbfact++]=e;}
  1202.       f1=f2;df1=df2;
  1203.     }
  1204.                                           /*printf("fin 1re etape\n");*/
  1205.   while(lgef(f1)>3);
  1206.   
  1207.                                     /*printf("fin de l'algorithme\n");*/
  1208.   nbf=nbfact;
  1209.   tetpil=avma;
  1210.   t[1]=gdiv(t[1],((GEN)t[1])[lgef(t[1])-1]);
  1211.   for(j=2;j<nbfact;j++)
  1212.     {
  1213.       if (expos[j]) t[j]=gdiv(t[j],((GEN)t[j])[lgef(t[j])-1]);
  1214.       for (k=1;k<j;k++)
  1215.     if (expos[k]&&polegal(t[j],t[k])) {expos[k]+=expos[j];expos[j]=0;nbf--;k=j;}
  1216.     }
  1217.   y=cgetg(3,19);
  1218.   u=cgetg(nbf,18);y[1]=(long)u;
  1219.   v=cgetg(nbf,18);y[2]=(long)v;
  1220.   for(j=1,k=0;j<nbfact;j++)
  1221.     {
  1222.       if (expos[j])
  1223.         {
  1224.           k++;
  1225.           u[k]=(long)t[j];
  1226.           v[k]=lstoi(expos[j]);
  1227.     }
  1228.     }
  1229.   return(gerepile(av,tetpil,y));
  1230. }
  1231.  
  1232.  
  1233. GEN stopoly(m,p,v) long m,p,v;
  1234. /* renvoie un polynome de la variable v
  1235. dont les coef sont les digits de m en base p */
  1236. {
  1237.   GEN y;long l=2,i,c[1000];
  1238.   do {c[l++]=m%p;m=m/p;} while(m);
  1239.   y=cgetg(l,10);for(i=2;i<l;i++) y[i]=lstoi(c[i]);
  1240.   y[1]=0x1000000+l;setvarn(y,v);
  1241.   return y;
  1242. }
  1243.  
  1244. void split(m,t,d,p,q)
  1245.      GEN *t,q;
  1246.      long p,d,m;
  1247.  
  1248. /*----------------------------------------------------- 
  1249.    Programme recursif :
  1250.   Entree:
  1251.    m entier arbitraire ( converti en un polynome w )
  1252.    p nb premier; q=(p^d-1)/2
  1253.    t[0] polynome de degre k*d prod de k fact de deg d.
  1254.   Sortie:
  1255.    t[0],t[1]...t[k-1] contiennent les k facteurs de g
  1256. ------------------------------------------------------*/
  1257.  
  1258. {  
  1259.   long j,l,v,av,av1,dv;
  1260.   GEN w,w0,wm,unmodp;
  1261.  
  1262.   if ((dv=lgef(*t)-3)==d) return;
  1263.   v=varn(*t);
  1264.   unmodp=gmodulcp(un,stoi(p));
  1265.   do
  1266.     {
  1267.       av=avma;
  1268.       if(p==2)
  1269.     {
  1270.       w=gmul(gpuigs(polx[v],m-1),unmodp);m+=2;
  1271.       for(w0=w,j=1;j<d;j++) w=gmod(gadd(w0,gmul(w,w)),*t);
  1272.     }
  1273.       else
  1274.     {
  1275.       w=gmul(stopoly(m++,p,v),unmodp);
  1276.       wm=gpui(gmodulcp(w,*t),q);w=gsub(wm[2],unmodp);
  1277.     }
  1278.       av1=avma;
  1279.       w=gerepile(av,av1,ggcd(*t,w));
  1280.       l=lgef(w)-3;
  1281.     }
  1282.   while((!l)||(l==dv));
  1283.   t[j=l/d]=gdiv(*t,w);*t=w;
  1284.   split(m,t+j,d,p,q);split(m,t,d,p,q);
  1285. }
  1286.  
  1287. GEN decpol(x,klim)
  1288.   GEN x;
  1289.   long klim;
  1290.   
  1291. /* a modifier. Ecrit principalement pour le programme trace.c */
  1292. /* aucune verification */
  1293. /* klim=0 habituellement, sauf si l'on ne veut chercher que les facteurs de degre <= klim */
  1294.  
  1295. {
  1296.   short int pos[200];
  1297.   long av=avma,av1,tete,lx,k,kin,i,j,i1,i2,fl,d,nbfact;
  1298.   GEN res,p1,p2;
  1299.   
  1300.   kin=1;res=cgetg(lx=lgef(x)-2,17);nbfact=0;
  1301.   p1=roots(x,4);d=lg(p1)-1;if(!klim) klim=d;
  1302.   do
  1303.   {
  1304.     fl=1;
  1305.     for(k=kin;((k+k)<=d)&&fl&&(k<=klim);k++)
  1306.     {
  1307.       for(i=0;i<=k;i++) pos[i]=i;
  1308.       do
  1309.       {
  1310.         av1=avma;p2=gzero;j=0;
  1311.     for(i1=1;i1<=k;i1++) p2=gadd(p2,p1[pos[i1]]);
  1312.         if((gexpo(gimag(p2))<-20)&&(gexpo(gsub(p2,ground(p2)))<-20))
  1313.         {
  1314.           p2=gun;
  1315.       for(i1=1;i1<=k;i1++) p2=gmul(p2,gsub(polx[0],p1[pos[i1]]));p2=ground(p2);
  1316.           if((gcmp0(gimag(p2)))&&(gcmp0(gmod(x,p2)))) 
  1317.           {
  1318.             res[++nbfact]=(long)p2;x=gdiv(x,p2);fl=0;kin=k;p2=cgetg(d-k+1,18);
  1319.             i1=1;i2=1;for(i=1;i<=d;i++)
  1320.             {
  1321.               if((i1<=k)&&(i==pos[i1])) i1++;
  1322.               else p2[i2++]=p1[i];
  1323.             }
  1324.             p1=p2;d-=k;
  1325.           }
  1326.         }
  1327.         if(fl)
  1328.         {
  1329.           avma=av1;pos[k]++;
  1330.           while(pos[k-j]>(d-j))
  1331.           {
  1332.             j++;pos[k-j]++;
  1333.           }
  1334.           for(i=k-j+1;i<=k;i++) pos[i]=i+pos[k-j]-k+j;
  1335.         }
  1336.       }
  1337.       while((j<k)&&fl);
  1338.     }
  1339.   }
  1340.   while(((((k+k)<=d)&&(k<=klim))||(!fl))&&(lgef(x)>3));
  1341.   if(lgef(x)>3) res[++nbfact]=(long)x;
  1342.   setlg(res,nbfact+1);for(j=nbfact+1;j<lx;j++) res[j]=zero; /*necessaire pour gerepile */
  1343.   tete=avma;return gerepile(av,tete,greal(res));
  1344. }
  1345.  
  1346.  
  1347. GEN factpol2(x,klim)
  1348.   GEN x;
  1349.   long klim;
  1350.  
  1351. /* klim=0 habituellement, sauf si l'on ne veut chercher que les facteurs de degre <= klim */
  1352.  
  1353. {
  1354.   long av=avma,av2,vv,k,i,j,i1,f,nbfac;
  1355.   GEN res,p1,p2,y,d,fa[30],a,ap,t,v,w;
  1356.   
  1357.   if((typ(x)!=10)||(!signe(x))) err(factpoler1);
  1358.   y=cgetg(3,19);if(lgef(x)==3) {y[1]=lgetg(1,18);y[2]=lgetg(1,18);return y;}
  1359.   if(lgef(x)==4)
  1360.     {
  1361.       p1=cgetg(2,18);y[1]=(long)p1;p1[1]=lcopy(x);
  1362.       p1=cgetg(2,18);y[2]=(long)p1;p1[1]=un;return y;
  1363.     }
  1364.   d=content(x);vv=varn(x);a=gdiv(x,d);ap=deriv(a,vv);t=ggcd(a,ap);v=gdiv(a,t);
  1365.   w=gdiv(ap,t);j=0;f=1;nbfac=0;
  1366.   while(f)
  1367.     {
  1368.       j++;w=gsub(w,deriv(v,vv));f=signe(w);
  1369.       if(f) {res=ggcd(v,w);v=gdiv(v,res);w=gdiv(w,res);}
  1370.       else res=v;
  1371.       fa[j]=(lgef(res)>3) ? decpol(res,klim) : cgetg(1,18);
  1372.       nbfac+=(lg(fa[j])-1);
  1373.     }
  1374.   av2=avma;y=cgetg(3,19);p1=cgetg(nbfac+1,18);y[1]=(long)p1;
  1375.   p2=cgetg(nbfac+1,18);y[2]=(long)p2;
  1376.   for(i=1,k=0;i<=j;i++)
  1377.     for(i1=1;i1<lg(fa[i]);i1++)
  1378.       {
  1379.     p1[++k]=lcopy(fa[i][i1]);p2[k]=lstoi(i);
  1380.       }
  1381.   return gerepile(av,av2,y);
  1382. }
  1383.  
  1384. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1385. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1386. /*                                                                 */
  1387. /*                Recherche de racines  p-adiques                  */
  1388. /*                                                                 */
  1389. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1390. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1391.  
  1392. /* a etant un p-adique, retourne le vecteur des racines p-adiques de f
  1393. congrues a a modulo p dans le cas ou on suppose f(a) congru a 0 modulo p
  1394. (ou a 4 si p=2). */
  1395.  
  1396. #define gmaxval(x,y) (gcmp0(x)?BIGINT:ggval(x,y))
  1397.  
  1398. GEN apprgen(f,a)
  1399.      GEN f,a;
  1400.  
  1401. {
  1402.   GEN fp,fg,p1,p,pro,idiot,idiot2,u,ip,quatre;
  1403.   long av=avma,tetpil,v,vv,ps,i,j,k,lu,n,fl2;
  1404.  
  1405.   if((typ(f)!=10)||(typ(a)!=7)||gcmp0(f)) err(rootper1);
  1406.   v=varn(f);fp=deriv(f,v);fg=ggcd(f,fp);
  1407.   if(lgef(fg)>3) {f=gdiv(f,fg);fp=deriv(f,v);}
  1408.   p=(GEN)a[2];p1=poleval(f,a);if((vv=gmaxval(p1,p))<=0) err(rootper2);
  1409.   fl2=gegal(p,gdeux);if((vv==1)&&fl2) err(rootper2);
  1410.   if(fl2) quatre=stoi(4);vv=gmaxval(poleval(fp,a),p);
  1411.   if(!vv)
  1412.     {
  1413.       while(!gcmp0(p1)) {a=gsub(a,gdiv(p1,poleval(fp,a)));p1=poleval(f,a);}
  1414.       tetpil=avma;pro=cgetg(2,17);pro[1]=lcopy(a);
  1415.       return gerepile(av,tetpil,pro);
  1416.     }
  1417.   n=lgef(f)-3;pro=cgetg(n+1,17);j=0;
  1418.   p1=poleval(f,gadd(a,gmul(fl2?quatre:p,polx[v])));
  1419.   if(!gcmp0(p1)) p1=gdiv(p1,gpuigs(p,ggval(p1,p)));
  1420.   if(gcmpgs(p,VERYBIGINT)>0) err(impl,"apprgen for p>=2^31");
  1421.   ps=fl2?4:itos(p);idiot=gsub(a,a);idiot2=fl2 ? ggrando(p,2) : ggrando(p,1);
  1422.   for(i=0;i<ps;i++)
  1423.     {
  1424.       ip=stoi(i);
  1425.       if(gcmp0(poleval(p1,gadd(ip,idiot2))))
  1426.     {
  1427.       u=apprgen(p1,gadd(idiot,ip));
  1428.       lu=lg(u);
  1429.       for(k=1;k<lu;k++) {j++;pro[j]=ladd(a,gmul(fl2?quatre:p,u[k]));}
  1430.     }
  1431.     }
  1432.   tetpil=avma;p1=cgetg(j+1,17);for(i=1;i<=j;i++) p1[i]=lcopy(pro[i]);
  1433.   return gerepile(av,tetpil,p1);
  1434. }
  1435.  
  1436. /* Retourne le vecteur des racines p-adiques de f en precision r */
  1437.  
  1438. GEN rootpadic(f,p,r)
  1439.      GEN f,p;
  1440.      long r;
  1441. {
  1442.   GEN y,fp,fg,pro,yi,p1,rac;
  1443.   long v,lx,i,j,k,n,av=avma,tetpil,fl2;
  1444.  
  1445.   if((typ(f)!=10)||gcmp0(f)) err(rootper1);
  1446.   if(r<=0) err(rootper4);
  1447.   v=varn(f);fp=deriv(f,v);fg=ggcd(f,fp);
  1448.   if(lgef(fg)>3) {f=gdiv(f,fg);fp=deriv(f,v);}
  1449.   fl2=gegal(p,gdeux);rac=(fl2&&(r>=2))? rootmod(f,stoi(4)) : rootmod(f,p);
  1450.   lx=lg(rac);
  1451.   if(r==1)
  1452.     {
  1453.       tetpil=avma;y=cgetg(lx,18);
  1454.       for(i=1;i<lx;i++)
  1455.     {
  1456.       p1=(GEN)rac[i];yi=cgetg(5,7);yi[2]=lclone(p);yi[3]=lcopy(p);
  1457.       yi[4]=lcopy(p1[2]);yi[1]=0x18000;y[i]=(long)yi;
  1458.     }
  1459.       return gerepile(av,tetpil,y);
  1460.     }
  1461.   n=lgef(f)-3;pro=cgetg(n+1,18);j=0;
  1462.   for(i=1;i<lx;i++)
  1463.     {
  1464.       p1=(GEN)rac[i];yi=cgetg(5,7);yi[2]=lclone(p);
  1465.       if(signe(p1[2]))
  1466.     {
  1467.       if((!fl2)||mpodd(p1[2]))
  1468.         {
  1469.           yi[3]=lpuigs(p,r);yi[4]=p1[2];yi[1]=0x8000+(r<<16);
  1470.         }
  1471.       else
  1472.         {
  1473.           yi[3]=lpuigs(p,r);yi[4]=un;yi[1]=0x8001+(r<<16);
  1474.         }
  1475.     }
  1476.       else {yi[3]=un;yi[4]=zero;yi[1]=0x8000+r;}
  1477.       p1=apprgen(f,yi);for(k=1;k<lg(p1);k++) pro[++j]=p1[k];
  1478.     }
  1479.   tetpil=avma;p1=cgetg(j+1,18);for(i=1;i<=j;i++) p1[i]=lcopy(pro[i]);
  1480.   return gerepile(av,tetpil,p1);
  1481. }  
  1482.  
  1483. /* a appartenant a une extension finie de Q_p, retourne le vecteur des racines
  1484. de f congrues a a modulo p dans le cas ou on suppose f(a) congru a 0 modulo p
  1485. (ou a 4 si p=2). */
  1486.  
  1487. GEN apprgen9(f,a)
  1488.      GEN f,a;
  1489.  
  1490. {
  1491.   GEN fp,fg,p1,p,pro,idiot,idiot2,u,ip,alpha,t,vecg,quatre;
  1492.   long av=avma,tetpil,v,vv,ps,i,j,k,lu,n,precpadique,d,va,flfin,fl2;
  1493.  
  1494.   if((typ(f)!=10)||gcmp0(f)) err(rootper1);
  1495.   if(typ(a)==7) return apprgen(f,a);
  1496.   if((typ(a)!=9)||(typ(a[2])!=10)) err(rootper1);
  1497.   v=varn(f);fp=deriv(f,v);fg=ggcd(f,fp);
  1498.   if(lgef(fg)>3) {f=gdiv(f,fg);fp=deriv(f,v);}
  1499.   alpha=(GEN)a[2];t=(GEN)a[1];precpadique=BIGINT;va=varn(t);
  1500.   for(i=2;i<lgef(alpha);i++)
  1501.     {
  1502.       pro=(GEN)alpha[i];
  1503.       if(typ(pro)==7) 
  1504.     {
  1505.       precpadique=min(precpadique,(signe(pro[4])?valp(pro)+precp(pro):valp(pro)));
  1506.       p=(GEN)pro[2];
  1507.     }
  1508.     }
  1509.   d=lgef(t)-3;
  1510.   for(i=2;i<=d+2;i++)
  1511.     {
  1512.       pro=(GEN)t[i];
  1513.       if(typ(pro)==7) 
  1514.     {
  1515.       precpadique=min(precpadique,(signe(pro[4])?valp(pro)+precp(pro):valp(pro)));
  1516.       p=(GEN)pro[2];
  1517.     }
  1518.     }
  1519.   if(precpadique==BIGINT) err(rootper1);
  1520.   p1=poleval(f,a);if((vv=gmaxval(lift(p1),p))<=0) err(rootper2);
  1521.   fl2=gegal(p,gdeux);if((vv==1)&&fl2) err(rootper2);
  1522.   if(fl2) quatre=stoi(4);vv=gmaxval(lift(poleval(fp,a)),p);
  1523.   if(!vv)
  1524.     {
  1525.       while(!gcmp0(p1)) {a=gsub(a,gdiv(p1,poleval(fp,a)));p1=poleval(f,a);}
  1526.       tetpil=avma;pro=cgetg(2,18);pro[1]=lcopy(a);
  1527.       return gerepile(av,tetpil,pro);
  1528.     }
  1529.   n=lgef(f)-3;pro=cgetg(n+1,18);j=0;
  1530.   p1=poleval(f,gadd(a,gmul(fl2?quatre:p,polx[v])));
  1531.   if(!gcmp0(p1)) p1=gdiv(p1,gpuigs(p,ggval(p1,p)));
  1532.   if(gcmpgs(p,VERYBIGINT)>0) err(impl,"apprgen9 for p>=2^31");
  1533.   ps=fl2?4:itos(p);idiot=gmodulcp(ggrando(p,precpadique),t);
  1534.   idiot2=gmodulcp(fl2 ? ggrando(p,2) : ggrando(p,1),t);
  1535.   vecg=cgetg(d+1,18);for(i=1;i<=d;i++) vecg[i]=zero;
  1536.   flfin=1;
  1537.   while(flfin)
  1538.     {
  1539.       ip=gmodulcp(gtopoly(vecg,va),t);
  1540.       if(gcmp0(poleval(p1,gadd(ip,idiot2))))
  1541.     {
  1542.       u=apprgen9(p1,gadd(ip,idiot));
  1543.       lu=lg(u);
  1544.       for(k=1;k<lu;k++) {j++;pro[j]=ladd(a,gmul(fl2?quatre:p,u[k]));}
  1545.     }
  1546.       i=d;while(i&&(!cmpis(vecg[i],ps-1))) i--;
  1547.       if(i) vecg[i]=laddsi(1,vecg[i]);
  1548.       else flfin=0;
  1549.     }
  1550.   tetpil=avma;p1=cgetg(j+1,18);for(i=1;i<=j;i++) p1[i]=lcopy(pro[i]);
  1551.   return gerepile(av,tetpil,p1);
  1552. }
  1553.  
  1554. GEN padicff(f,p,r)
  1555.      GEN f,p;
  1556.      long r;
  1557.  
  1558. /* interne donc aucune verification */
  1559. /* En entree, res est un polynome a coeffs dans Z_p sans facteur carre */
  1560.  
  1561. {
  1562.   long av=avma,tetpil,lx,n,i,j,k,v,fl2;
  1563.   GEN rac,pro,p1,p2,y,yi,quatre;
  1564.  
  1565.   fl2=gegal(p,gdeux);
  1566.   if(fl2&&(r>=2)) err(impl,"padicfactor for p=2");
  1567.   rac=factmod(f,p);lx=lg(rac[1]);
  1568.   if(r==1)
  1569.     {
  1570.       tetpil=avma;y=cgetg(lx,18);p2=(GEN)rac[1];
  1571.       for(i=1;i<lx;i++) y[i]=(long)gcvtop(p2[i],p,1);
  1572.       return gerepile(av,tetpil,y);
  1573.     }
  1574.   if(fl2) quatre=stoi(4);
  1575.   n=lgef(f)-3;pro=cgetg(n+1,18);j=0;v=varn(f);p2=lift(rac[1]);
  1576.   for(i=1;i<lx;i++)
  1577.     {
  1578.       p1=(GEN)p2[i];
  1579.       if(lgef(p1)==4)
  1580.     {
  1581.       p1=gcmp0(p1[2])?(GEN)p1[2]:gsub(fl2?quatre:p,p1[2]);
  1582.       yi=cgetg(5,7);yi[2]=lclone(p);
  1583.       if(signe(p1))
  1584.         {
  1585.           if((!fl2)||mpodd(p1))
  1586.         {
  1587.           yi[3]=lpuigs(p,r);yi[4]=(long)p1;yi[1]=0x8000+(r<<16);
  1588.         }
  1589.           else
  1590.         {
  1591.           yi[3]=lpuigs(p,r);yi[4]=un;yi[1]=0x8001+(r<<16);
  1592.         }
  1593.         }
  1594.       else {yi[3]=un;yi[4]=zero;yi[1]=0x8000+r;}
  1595.       p1=apprgen(f,yi);
  1596.       for(k=1;k<lg(p1);k++) pro[++j]=lsub(polx[v],p1[k]);
  1597.     }
  1598.       else
  1599.     {
  1600.       yi=gmodulcp(gcvtop(polx[v],p,r),gcvtop(p2[i],p,r));
  1601.       p1=apprgen9(f,yi);
  1602.       for(k=1;k<lg(p1);k++) pro[++j]=(long)caract(p1[k],v);
  1603.     }
  1604.     }
  1605.   tetpil=avma;p1=cgetg(j+1,18);for(i=1;i<=j;i++) p1[i]=lcopy(pro[i]);
  1606.   return gerepile(av,tetpil,p1);
  1607. }  
  1608.  
  1609. GEN factorpadic(x,p,r)
  1610.   GEN x,p;
  1611.   long r;
  1612.  
  1613. {
  1614.   long av=avma,av2,vv,k,i,j,i1,f,nbfac;
  1615.   GEN res,p1,p2,y,d,fa[100],a,ap,t,v,w;
  1616.   
  1617.   if((typ(x)!=10)||(!signe(x))) err(rootper1);
  1618.   if(r<=0) err(rootper4);
  1619.   y=cgetg(3,19);if(lgef(x)==3) {y[1]=lgetg(1,18);y[2]=lgetg(1,18);return y;}
  1620.   if(lgef(x)==4)
  1621.     {
  1622.       p1=cgetg(2,18);y[1]=(long)p1;p1[1]=lcopy(x);
  1623.       p1=cgetg(2,18);y[2]=(long)p1;p1[1]=un;return y;
  1624.     }
  1625.   d=content(x);vv=varn(x);a=gdiv(x,d);ap=deriv(a,vv);t=ggcd(a,ap);v=gdiv(a,t);
  1626.   w=gdiv(ap,t);j=0;f=1;nbfac=0;
  1627.   while(f)
  1628.     {
  1629.       j++;w=gsub(w,deriv(v,vv));f=signe(w);
  1630.       if(f)
  1631.     {
  1632.       res=ggcd(v,w);v=gdiv(v,res);w=gdiv(w,res);
  1633.     }
  1634.       else res=v;
  1635.       fa[j]=(lgef(res)>3) ? padicff(res,p,r) : cgetg(1,18);
  1636.       nbfac+=(lg(fa[j])-1);
  1637.     }
  1638.   av2=avma;y=cgetg(3,19);p1=cgetg(nbfac+1,18);y[1]=(long)p1;
  1639.   p2=cgetg(nbfac+1,18);y[2]=(long)p2;
  1640.   for(i=1,k=0;i<=j;i++)
  1641.     for(i1=1;i1<lg(fa[i]);i1++)
  1642.       {
  1643.     p1[++k]=lcopy(fa[i][i1]);p2[k]=lstoi(i);
  1644.       }
  1645.   return gerepile(av,av2,y);
  1646. }
  1647.  
  1648.  
  1649. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1650. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1651. /*                                                                 */
  1652. /*                     FACTORISATION DANS F_q                      */
  1653. /*                                                                 */
  1654. /*                                                                 */
  1655. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1656. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1657.  
  1658. GEN     factmod9(f,p,a)
  1659.      
  1660.      GEN     f,p,a;
  1661.      
  1662. {
  1663.   long  i,j,k,d,e,vf,va,expos[100],psim,nbfact,nbf,av,pk,tetpil,calc,qqs;
  1664.   GEN   y,t[100],f1,f2,f3,df1,df2,g,g1,g2;
  1665.   GEN   xmod,u,v,pd,q,qq,unfp,unfq;
  1666.   
  1667.   if((typ(a)!=10)||(typ(f)!=10)||gcmp0(f)||gcmp0(a)) err(factmoder);
  1668.   if(lgef(f)==3)
  1669.     {y=cgetg(3,19);y[1]=lgetg(1,18);y[2]=lgetg(1,18);return y;}
  1670.   vf=varn(f);va=varn(a);av=avma;unfp=gmodulcp(un,p);
  1671.   a=gmul(unfp,a);unfq=gmodulcp(gmul(unfp,polun[va]),a);
  1672.   f=gmul(unfq,f);qq=gpuigs(p,lgef(a)-3);qqs=itos(qq);
  1673.   
  1674.   /*  1ere etape : trouver les facteurs square-free produit
  1675.       des facteurs premiers de meme multiplicite    */
  1676.  
  1677.   f1=f;e=0;nbfact=1;pk=1;df1=deriv(f1,vf);calc=0;
  1678.  
  1679.   do
  1680.     {
  1681. /*printf("debut 1re etape\n");*/
  1682.       e+=pk;
  1683.       while (gcmp0(df1))
  1684.         {
  1685.           pk *= (psim=itos(p));e=pk;
  1686.           j=(lgef(f1)-3)/psim+3;f2=cgetg(j,10);
  1687.           f2[1]=0x1000000+j;setvarn(f2,vf);
  1688.           for(i=2;i<j;i++)
  1689.             f2[i]=lcopy(f1[psim*(i-2)+2]);
  1690.           f1=f2;df1=deriv(f1,vf);calc=0;
  1691.     }
  1692.       if(calc) f2=f3;else f2=ggcd(f1,df1);calc=1;
  1693.       if (lgef(f2)<4) u=f1;
  1694.       else
  1695.         {
  1696.           g1=gdiv(f1,f2);
  1697.           if (gcmp0(df2=deriv(f2,vf)))
  1698.             {
  1699.               u=g1;f3=f2;
  1700.         }
  1701.           else
  1702.             {
  1703.               f3=ggcd(f2,df2);
  1704.               if (lgef(f3)<4) u=gdiv(g1,f2);
  1705.               else
  1706.                 {
  1707.                   g2=gdiv(f2,f3);u=gdiv(g1,g2);
  1708.          }
  1709.         }
  1710.     }
  1711.       
  1712.       /*  2eme etape : 
  1713.           Ici u est un polynome square-free
  1714.           (produit des facteurs premiers de meme multiplicite  e :
  1715.           trouver les facteurs (square-free) produit des facteurs de meme
  1716.       degre d */
  1717.       
  1718.       d=0;pd=gun;
  1719.       xmod=gmodulcp(polx[vf],u);v=xmod;
  1720.       while(d<(lgef(u)-3)>>1)
  1721.         {
  1722.                                        /*printf("debut 2me etape\n");*/
  1723.           d++;
  1724.           pd=mulii(pd,qq);
  1725.           q=shifti(subis(pd,1),-1);
  1726.           v=gpui(v,qq);
  1727.           g=ggcd((gsub(v,xmod))[2],u);
  1728.       
  1729.           if (lgef(g)>3)
  1730.             {
  1731.                                       /*printf("debut 3me etape\n");*/
  1732.           
  1733. /*  3eme etape :
  1734. Ici g est produit de pol irreductibles ayant tous le meme degre d;*/
  1735.  
  1736.           t[nbfact]=g;j=nbfact+(lgef(g)-3)/d;
  1737.           split9(qqs,t+nbfact,d,p[2],q,unfq,qqs,a);
  1738. /* le premier parametre est un entier variable m qui sera converti en un
  1739.    polynome w dont les coeff sont ses digits en base p (initialement m = p --> X)
  1740.    pour faire pgcd de g avec w^(qq^d-1)/2 jusqu'a casser. */
  1741.           for(;nbfact<j;expos[nbfact++]=e);
  1742.           u=gdiv(u,g);v=gmodulcp(v[2],u);
  1743.         }
  1744.                                           /*printf("fin 3me etape\n");*/
  1745.  
  1746.     }                                 /*printf("fin 2me etape\n");*/
  1747.       if (lgef(u)>3) {t[nbfact]=u;expos[nbfact++]=e;}
  1748.       f1=f2;df1=df2;
  1749.     }
  1750.                                           /*printf("fin 1re etape\n");*/
  1751.   while(lgef(f1)>3);
  1752.   
  1753.                                     /*printf("fin de l'algorithme\n");*/
  1754.   nbf=nbfact;
  1755.   tetpil=avma;
  1756.   t[1]=gdiv(t[1],((GEN)t[1])[lgef(t[1])-1]);
  1757.   for(j=2;j<nbfact;j++)
  1758.     {
  1759.       if (expos[j]) t[j]=gdiv(t[j],((GEN)t[j])[lgef(t[j])-1]);
  1760.       for (k=1;k<j;k++)
  1761.     if (expos[k]&&polegal(t[j],t[k])) {expos[k]+=expos[j];expos[j]=0;nbf--;k=j;}
  1762.     }
  1763.   y=cgetg(3,19);
  1764.   u=cgetg(nbf,18);y[1]=(long)u;
  1765.   v=cgetg(nbf,18);y[2]=(long)v;
  1766.   for(j=1,k=0;j<nbfact;j++)
  1767.     {
  1768.       if (expos[j])
  1769.         {
  1770.           k++;
  1771.           u[k]=(long)t[j];
  1772.           v[k]=lstoi(expos[j]);
  1773.     }
  1774.     }
  1775.   return(gerepile(av,tetpil,y));
  1776. }
  1777.  
  1778. GEN stopoly9(m,p,qqs,v,a)
  1779.      long p,v,m,qqs;
  1780.      GEN a;
  1781.  
  1782. /* renvoie un polynome de la variable v
  1783. dont les coef sont les digits de m en base qq */
  1784.  
  1785. {
  1786.   GEN y,p2;
  1787.   long p1,l=2,l1,i,j,va=varn(a),c[1000],d[1000];
  1788.  
  1789.   do {c[l++]=m%qqs;m/=qqs;} while(m);
  1790.   y=cgetg(l,10);y[1]=0x1000000+l;setvarn(y,v);
  1791.   for(i=2;i<l;i++) 
  1792.     {
  1793.       p1=c[i];l1=2;
  1794.       do {d[l1++]=p1%p;p1/=p;} while(p1);
  1795.       p2=cgetg(l1,10);p2[1]=0x1000000+l1;setvarn(p2,va);
  1796.       for(j=2;j<l1;j++) p2[j]=lstoi(d[j]);
  1797.       y[i]=lmodulcp(p2,a);
  1798.     }
  1799.   return y;
  1800. }
  1801.  
  1802. GEN stopoly92(d1,v,a,ptres)
  1803.      long v,d1;
  1804.      GEN a,*ptres;
  1805.  
  1806. /* renvoie un polynome aleatoire de la variable v
  1807. de degre inferieur ou egal a 2*d1-1 */
  1808.  
  1809. {
  1810.   GEN y,p2;
  1811.   long p1,l1,i,j,d2,l,va=varn(a),c[1000],d[1000],k=lgef(a)-3;
  1812.  
  1813. /* qqs=2^k */
  1814.  
  1815.   c[1]=1;d2=d1+d1+1;
  1816.   do
  1817.     {
  1818.       for(l=2;l<=d2;l++) c[l]=(rand())>>(31-k);
  1819.       l=d2;while(!c[l]) l--;
  1820.     }
  1821.   while(l<=2);
  1822.   l++;
  1823.   y=cgetg(l,10);y[1]=0x1000000+l;setvarn(y,v);
  1824.   for(i=2;i<l;i++) 
  1825.     {
  1826.       p1=c[i];l1=2;
  1827.       do {d[l1++]=p1&1;p1>>=1;} while(p1);
  1828.       p2=cgetg(l1,10);p2[1]=0x1000000+l1;setvarn(p2,va);
  1829.       for(j=2;j<l1;j++) p2[j]=lstoi(d[j]);
  1830.       y[i]=lmodulcp(p2,a);
  1831.     }
  1832.   p1=(rand())>>(31-k);l1=2;
  1833.   do {d[l1++]=p1&1;p1>>=1;} while(p1);
  1834.   p2=cgetg(l1,10);p2[1]=0x1000000+l1;setvarn(p2,va);
  1835.   for(j=2;j<l1;j++) p2[j]=lstoi(d[j]);
  1836.   *ptres=gmodulcp(p2,a);
  1837.   return y;
  1838. }
  1839.  
  1840. void split9(m,t,d,p,q,unfq,qqs,a)
  1841.      GEN *t,q,unfq,a;
  1842.      long m,d,p,qqs;
  1843.  
  1844. /*----------------------------------------------------- 
  1845.    Programme recursif :
  1846.   Entree:
  1847.    m entier arbitraire ( converti en un polynome w )
  1848.    p nb premier; q=(qq^d-1)/2
  1849.    t[0] polynome de degre k*d prod de k fact de deg d.
  1850.   Sortie:
  1851.    t[0],t[1]...t[k-1] contiennent les k facteurs de g
  1852. ------------------------------------------------------*/
  1853.  
  1854. {  
  1855.   long j,l,v,av,av1,dv;
  1856.   GEN w,w0,wm,res;
  1857.  
  1858.   if ((dv=lgef(*t)-3)==d) return;
  1859.   v=varn(*t);
  1860.   do
  1861.     {
  1862.       av=avma;
  1863.       if(p==2)
  1864.     {
  1865.       w=gmul(stopoly92(d,v,a,&res),unfq);
  1866.       for(w0=w,j=1;j<d;j++) w=gmod(gadd(w0,gpuigs(w,qqs)),*t);
  1867.       w=gsub(w,res);
  1868.     }
  1869.       else
  1870.     {
  1871.       w=gmul(stopoly9(m,p,qqs,v,a),unfq);m++;
  1872.       wm=gpui(gmodulcp(w,*t),q);w=gsub(wm[2],unfq);
  1873.     }
  1874.       av1=avma;
  1875.       w=gerepile(av,av1,ggcd(*t,w));
  1876.       l=lgef(w)-3;
  1877.     }
  1878.   while((!l)||(l==dv));
  1879.   t[j=l/d]=gdiv(*t,w);*t=w;
  1880.   split9(m,t+j,d,p,q,unfq,qqs,a);split9(m,t,d,p,q,unfq,qqs,a);
  1881. }
  1882.  
  1883. GEN randompoly9(d1,p,v,a)
  1884.      long v,d1;
  1885.      GEN a,p;
  1886.  
  1887. /* renvoie un polynome aleatoire de la variable v
  1888. de degre inferieur ou egal a d1 */
  1889.  
  1890. {
  1891.   GEN y,p2;
  1892.   long p1,ps,qqs,l1,i,j,l,va=varn(a),c[1000],d[1000],k=lgef(a)-3;
  1893.  
  1894.   c[1]=1;qqs=itos(gpuigs(p,k));ps=itos(p);
  1895.   do
  1896.     {
  1897.       for(l=2;l<=d1+2;l++) c[l]=(rand()*(((double)qqs)/C31));
  1898.       l=d1+2;while(!c[l]) l--;
  1899.     }
  1900.   while(l<=2);
  1901.   l++;
  1902.   y=cgetg(l,10);y[1]=0x1000000+l;setvarn(y,v);
  1903.   for(i=2;i<l;i++) 
  1904.     {
  1905.       p1=c[i];l1=2;
  1906.       do {d[l1++]=p1%ps;p1/=ps;} while(p1);
  1907.       p2=cgetg(l1,10);p2[1]=0x1000000+l1;setvarn(p2,va);
  1908.       for(j=2;j<l1;j++) p2[j]=lstoi(d[j]);
  1909.       y[i]=lmodulcp(p2,a);
  1910.     }
  1911.   return y;
  1912. }
  1913.  
  1914.  
  1915.