home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 2 / RISC_DISC_2.iso / pd_share / program / code / pari / c / polarit1 < prev    next >
Encoding:
Text File  |  1991-11-28  |  51.5 KB  |  1,914 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. 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. 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.   gaffect(p1,z);avma=av;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;
  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.         *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(invmoder);
  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(gcmpgs(p2,0)>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. GEN stopoly(m,p,v) long m,p,v;
  1233. /* renvoie un polynome de la variable v
  1234. dont les coef sont les digits de m en base p */
  1235. {
  1236.   GEN y;long l=2,i,c[1000];
  1237.   do {c[l++]=m%p;m=m/p;} while(m);
  1238.   y=cgetg(l,10);for(i=2;i<l;i++) y[i]=lstoi(c[i]);
  1239.   y[1]=0x1000000+l;setvarn(y,v);
  1240.   return y;
  1241. }
  1242.  
  1243. split(m,t,d,p,q)
  1244.      GEN *t,q;
  1245.      long p,d,m;
  1246.  
  1247. /*----------------------------------------------------- 
  1248.    Programme recursif :
  1249.   Entree:
  1250.    m entier arbitraire ( converti en un polynome w )
  1251.    p nb premier; q=(p^d-1)/2
  1252.    t[0] polynome de degre k*d prod de k fact de deg d.
  1253.   Sortie:
  1254.    t[0],t[1]...t[k-1] contiennent les k facteurs de g
  1255. ------------------------------------------------------*/
  1256.  
  1257. {  
  1258.   long j,l,v,av,av1,dv;
  1259.   GEN w,w0,wm,unmodp;
  1260.  
  1261.   if ((dv=lgef(*t)-3)==d) return 0;
  1262.   v=varn(*t);
  1263.   unmodp=gmodulcp(un,stoi(p));
  1264.   do
  1265.     {
  1266.       av=avma;
  1267.       if(p==2)
  1268.     {
  1269.       w=gmul(gpuigs(polx[v],m-1),unmodp);m+=2;
  1270.       for(w0=w,j=1;j<d;j++) w=gmod(gadd(w0,gmul(w,w)),*t);
  1271.     }
  1272.       else
  1273.     {
  1274.       w=gmul(stopoly(m++,p,v),unmodp);
  1275.       wm=gpui(gmodulcp(w,*t),q);w=gsub(wm[2],unmodp);
  1276.     }
  1277.       av1=avma;
  1278.       w=gerepile(av,av1,ggcd(*t,w));
  1279.       l=lgef(w)-3;
  1280.     }
  1281.   while((!l)||(l==dv));
  1282.   t[j=l/d]=gdiv(*t,w);*t=w;
  1283.   split(m,t+j,d,p,q);split(m,t,d,p,q);
  1284. }
  1285.  
  1286. GEN decpol(x,klim)
  1287.   GEN x;
  1288.   long klim;
  1289.   
  1290. /* a modifier. Ecrit principalement pour le programme trace.c */
  1291. /* aucune verification */
  1292. /* klim=0 habituellement, sauf si l'on ne veut chercher que les facteurs de degre <= klim */
  1293.  
  1294. {
  1295.   short int pos[200];
  1296.   long av=avma,av1,tete,lx,k,kin,i,j,i1,i2,fl,d,nbfact;
  1297.   GEN res,p1,p2;
  1298.   
  1299.   kin=1;res=cgetg(lx=lgef(x)-2,17);nbfact=0;
  1300.   p1=roots(x,4);d=lg(p1)-1;if(!klim) klim=d;
  1301.   do
  1302.   {
  1303.     fl=1;
  1304.     for(k=kin;((k+k)<=d)&&fl&&(k<=klim);k++)
  1305.     {
  1306.       for(i=0;i<=k;i++) pos[i]=i;
  1307.       do
  1308.       {
  1309.         av1=avma;p2=gzero;j=0;
  1310.     for(i1=1;i1<=k;i1++) p2=gadd(p2,p1[pos[i1]]);
  1311.         if((gexpo(gimag(p2))<-20)&&(gexpo(gsub(p2,ground(p2)))<-20))
  1312.         {
  1313.           p2=gun;
  1314.       for(i1=1;i1<=k;i1++) p2=gmul(p2,gsub(polx[0],p1[pos[i1]]));p2=ground(p2);
  1315.           if((gcmp0(gimag(p2)))&&(gcmp0(gmod(x,p2)))) 
  1316.           {
  1317.             res[++nbfact]=(long)p2;x=gdiv(x,p2);fl=0;kin=k;p2=cgetg(d-k+1,18);
  1318.             i1=1;i2=1;for(i=1;i<=d;i++)
  1319.             {
  1320.               if((i1<=k)&&(i==pos[i1])) i1++;
  1321.               else p2[i2++]=p1[i];
  1322.             }
  1323.             p1=p2;d-=k;
  1324.           }
  1325.         }
  1326.         if(fl)
  1327.         {
  1328.           avma=av1;pos[k]++;
  1329.           while(pos[k-j]>(d-j))
  1330.           {
  1331.             j++;pos[k-j]++;
  1332.           }
  1333.           for(i=k-j+1;i<=k;i++) pos[i]=i+pos[k-j]-k+j;
  1334.         }
  1335.       }
  1336.       while((j<k)&&fl);
  1337.     }
  1338.   }
  1339.   while(((((k+k)<=d)&&(k<=klim))||(!fl))&&(lgef(x)>3));
  1340.   if(lgef(x)>3) res[++nbfact]=(long)x;
  1341.   setlg(res,nbfact+1);for(j=nbfact+1;j<lx;j++) res[j]=zero; /*necessaire pour gerepile */
  1342.   tete=avma;return gerepile(av,tete,greal(res));
  1343. }
  1344.  
  1345.  
  1346. GEN factpol2(x,klim)
  1347.   GEN x;
  1348.   long klim;
  1349.  
  1350. /* klim=0 habituellement, sauf si l'on ne veut chercher que les facteurs de degre <= klim */
  1351.  
  1352. {
  1353.   long av=avma,av2,vv,k,i,j,i1,f,nbfac;
  1354.   GEN res,p1,p2,y,d,fa[30],a,ap,t,v,w;
  1355.   
  1356.   if((typ(x)!=10)||(!signe(x))) err(factpoler1);
  1357.   y=cgetg(3,19);if(lgef(x)==3) {y[1]=lgetg(1,18);y[2]=lgetg(1,18);return y;}
  1358.   if(lgef(x)==4)
  1359.     {
  1360.       p1=cgetg(2,18);y[1]=(long)p1;p1[1]=lcopy(x);
  1361.       p1=cgetg(2,18);y[2]=(long)p1;p1[1]=un;return y;
  1362.     }
  1363.   d=content(x);vv=varn(x);a=gdiv(x,d);ap=deriv(a,vv);t=ggcd(a,ap);v=gdiv(a,t);
  1364.   w=gdiv(ap,t);j=0;f=1;nbfac=0;
  1365.   while(f)
  1366.     {
  1367.       j++;w=gsub(w,deriv(v,vv));f=signe(w);
  1368.       if(f) {res=ggcd(v,w);v=gdiv(v,res);w=gdiv(w,res);}
  1369.       else res=v;
  1370.       fa[j]=(lgef(res)>3) ? decpol(res,klim) : cgetg(1,18);
  1371.       nbfac+=(lg(fa[j])-1);
  1372.     }
  1373.   av2=avma;y=cgetg(3,19);p1=cgetg(nbfac+1,18);y[1]=(long)p1;
  1374.   p2=cgetg(nbfac+1,18);y[2]=(long)p2;
  1375.   for(i=1,k=0;i<=j;i++)
  1376.     for(i1=1;i1<lg(fa[i]);i1++)
  1377.       {
  1378.     p1[++k]=lcopy(fa[i][i1]);p2[k]=lstoi(i);
  1379.       }
  1380.   return gerepile(av,av2,y);
  1381. }
  1382.  
  1383. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1384. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1385. /*                                                                 */
  1386. /*                Recherche de racines  p-adiques                  */
  1387. /*                                                                 */
  1388. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1389. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1390.  
  1391. /* a etant un p-adique, retourne le vecteur des racines p-adiques de f
  1392. congrues a a modulo p dans le cas ou on suppose f(a) congru a 0 modulo p
  1393. (ou a 4 si p=2). */
  1394.  
  1395. #define gmaxval(x,y) (gcmp0(x)?BIGINT:ggval(x,y))
  1396.  
  1397. GEN apprgen(f,a)
  1398.      GEN f,a;
  1399.  
  1400. {
  1401.   GEN fp,fg,p1,p,pro,idiot,idiot2,u,ip,quatre;
  1402.   long av=avma,tetpil,v,vv,ps,i,j,k,lu,n,fl2;
  1403.  
  1404.   if((typ(f)!=10)||(typ(a)!=7)||gcmp0(f)) err(rootper1);
  1405.   v=varn(f);fp=deriv(f,v);fg=ggcd(f,fp);
  1406.   if(lgef(fg)>3) {f=gdiv(f,fg);fp=deriv(f,v);}
  1407.   p=(GEN)a[2];p1=poleval(f,a);if((vv=gmaxval(p1,p))<=0) err(rootper2);
  1408.   fl2=gegal(p,gdeux);if((vv==1)&&fl2) err(rootper2);
  1409.   if(fl2) quatre=stoi(4);vv=gmaxval(poleval(fp,a),p);
  1410.   if(!vv)
  1411.     {
  1412.       while(!gcmp0(p1)) {a=gsub(a,gdiv(p1,poleval(fp,a)));p1=poleval(f,a);}
  1413.       tetpil=avma;pro=cgetg(2,17);pro[1]=lcopy(a);
  1414.       return gerepile(av,tetpil,pro);
  1415.     }
  1416.   n=lgef(f)-3;pro=cgetg(n+1,17);j=0;
  1417.   p1=poleval(f,gadd(a,gmul(fl2?quatre:p,polx[v])));
  1418.   if(!gcmp0(p1)) p1=gdiv(p1,gpuigs(p,ggval(p1,p)));
  1419.   if(gcmpgs(p,VERYBIGINT)>0) err(impl,"apprgen for p>=2^31");
  1420.   ps=fl2?4:itos(p);idiot=gsub(a,a);idiot2=fl2 ? ggrando(p,2) : ggrando(p,1);
  1421.   for(i=0;i<ps;i++)
  1422.     {
  1423.       ip=stoi(i);
  1424.       if(gcmp0(poleval(p1,gadd(ip,idiot2))))
  1425.     {
  1426.       u=apprgen(p1,gadd(idiot,ip));
  1427.       lu=lg(u);
  1428.       for(k=1;k<lu;k++) {j++;pro[j]=ladd(a,gmul(fl2?quatre:p,u[k]));}
  1429.     }
  1430.     }
  1431.   tetpil=avma;p1=cgetg(j+1,17);for(i=1;i<=j;i++) p1[i]=lcopy(pro[i]);
  1432.   return gerepile(av,tetpil,p1);
  1433. }
  1434.  
  1435. /* Retourne le vecteur des racines p-adiques de f en precision r */
  1436.  
  1437. GEN rootpadic(f,p,r)
  1438.      GEN f,p;
  1439.      long r;
  1440. {
  1441.   GEN y,fp,fg,pro,yi,p1,rac;
  1442.   long v,lx,i,j,k,n,av=avma,tetpil,fl2;
  1443.  
  1444.   if((typ(f)!=10)||gcmp0(f)) err(rootper1);
  1445.   if(r<=0) err(rootper4);
  1446.   v=varn(f);fp=deriv(f,v);fg=ggcd(f,fp);
  1447.   if(lgef(fg)>3) {f=gdiv(f,fg);fp=deriv(f,v);}
  1448.   fl2=gegal(p,gdeux);rac=(fl2&&(r>=2))? rootmod(f,stoi(4)) : rootmod(f,p);
  1449.   lx=lg(rac);
  1450.   if(r==1)
  1451.     {
  1452.       tetpil=avma;y=cgetg(lx,18);
  1453.       for(i=1;i<lx;i++)
  1454.     {
  1455.       p1=(GEN)rac[i];yi=cgetg(5,7);yi[2]=lclone(p);yi[3]=lcopy(p);
  1456.       yi[4]=lcopy(p1[2]);yi[1]=0x18000;y[i]=(long)yi;
  1457.     }
  1458.       return gerepile(av,tetpil,y);
  1459.     }
  1460.   n=lgef(f)-3;pro=cgetg(n+1,18);j=0;
  1461.   for(i=1;i<lx;i++)
  1462.     {
  1463.       p1=(GEN)rac[i];yi=cgetg(5,7);yi[2]=lclone(p);
  1464.       if(signe(p1[2]))
  1465.     {
  1466.       if((!fl2)||mpodd(p1[2]))
  1467.         {
  1468.           yi[3]=lpuigs(p,r);yi[4]=p1[2];yi[1]=0x8000+(r<<16);
  1469.         }
  1470.       else
  1471.         {
  1472.           yi[3]=lpuigs(p,r);yi[4]=un;yi[1]=0x8001+(r<<16);
  1473.         }
  1474.     }
  1475.       else {yi[3]=un;yi[4]=zero;yi[1]=0x8000+r;}
  1476.       p1=apprgen(f,yi);for(k=1;k<lg(p1);k++) pro[++j]=p1[k];
  1477.     }
  1478.   tetpil=avma;p1=cgetg(j+1,18);for(i=1;i<=j;i++) p1[i]=lcopy(pro[i]);
  1479.   return gerepile(av,tetpil,p1);
  1480. }  
  1481.  
  1482. /* a appartenant a une extension finie de Q_p, retourne le vecteur des racines
  1483. de f congrues a a modulo p dans le cas ou on suppose f(a) congru a 0 modulo p
  1484. (ou a 4 si p=2). */
  1485.  
  1486. GEN apprgen9(f,a)
  1487.      GEN f,a;
  1488.  
  1489. {
  1490.   GEN fp,fg,p1,p,pro,idiot,idiot2,u,ip,alpha,t,vecg,quatre;
  1491.   long av=avma,tetpil,v,vv,ps,i,j,k,lu,n,precpadique,d,va,flfin,fl2;
  1492.  
  1493.   if((typ(f)!=10)||gcmp0(f)) err(rootper1);
  1494.   if(typ(a)==7) return apprgen(f,a);
  1495.   if((typ(a)!=9)||(typ(a[2])!=10)) err(rootper1);
  1496.   v=varn(f);fp=deriv(f,v);fg=ggcd(f,fp);
  1497.   if(lgef(fg)>3) {f=gdiv(f,fg);fp=deriv(f,v);}
  1498.   alpha=(GEN)a[2];t=(GEN)a[1];precpadique=BIGINT;va=varn(t);
  1499.   for(i=2;i<lgef(alpha);i++)
  1500.     {
  1501.       pro=(GEN)alpha[i];
  1502.       if(typ(pro)==7) 
  1503.     {
  1504.       precpadique=min(precpadique,(signe(pro[4])?valp(pro)+precp(pro):valp(pro)));
  1505.       p=(GEN)pro[2];
  1506.     }
  1507.     }
  1508.   d=lgef(t)-3;
  1509.   for(i=2;i<=d+2;i++)
  1510.     {
  1511.       pro=(GEN)t[i];
  1512.       if(typ(pro)==7) 
  1513.     {
  1514.       precpadique=min(precpadique,(signe(pro[4])?valp(pro)+precp(pro):valp(pro)));
  1515.       p=(GEN)pro[2];
  1516.     }
  1517.     }
  1518.   if(precpadique==BIGINT) err(rootper1);
  1519.   p1=poleval(f,a);if((vv=gmaxval(lift(p1),p))<=0) err(rootper2);
  1520.   fl2=gegal(p,gdeux);if((vv==1)&&fl2) err(rootper2);
  1521.   if(fl2) quatre=stoi(4);vv=gmaxval(lift(poleval(fp,a)),p);
  1522.   if(!vv)
  1523.     {
  1524.       while(!gcmp0(p1)) {a=gsub(a,gdiv(p1,poleval(fp,a)));p1=poleval(f,a);}
  1525.       tetpil=avma;pro=cgetg(2,18);pro[1]=lcopy(a);
  1526.       return gerepile(av,tetpil,pro);
  1527.     }
  1528.   n=lgef(f)-3;pro=cgetg(n+1,18);j=0;
  1529.   p1=poleval(f,gadd(a,gmul(fl2?quatre:p,polx[v])));
  1530.   if(!gcmp0(p1)) p1=gdiv(p1,gpuigs(p,ggval(p1,p)));
  1531.   if(gcmpgs(p,VERYBIGINT)>0) err(impl,"apprgen9 for p>=2^31");
  1532.   ps=fl2?4:itos(p);idiot=gmodulcp(ggrando(p,precpadique),t);
  1533.   idiot2=gmodulcp(fl2 ? ggrando(p,2) : ggrando(p,1),t);
  1534.   vecg=cgetg(d+1,18);for(i=1;i<=d;i++) vecg[i]=zero;
  1535.   flfin=1;
  1536.   while(flfin)
  1537.     {
  1538.       ip=gmodulcp(gtopoly(vecg,va),t);
  1539.       if(gcmp0(poleval(p1,gadd(ip,idiot2))))
  1540.     {
  1541.       u=apprgen9(p1,gadd(ip,idiot));
  1542.       lu=lg(u);
  1543.       for(k=1;k<lu;k++) {j++;pro[j]=ladd(a,gmul(fl2?quatre:p,u[k]));}
  1544.     }
  1545.       i=d;while(i&&(!cmpis(vecg[i],ps-1))) i--;
  1546.       if(i) vecg[i]=laddsi(1,vecg[i]);
  1547.       else flfin=0;
  1548.     }
  1549.   tetpil=avma;p1=cgetg(j+1,18);for(i=1;i<=j;i++) p1[i]=lcopy(pro[i]);
  1550.   return gerepile(av,tetpil,p1);
  1551. }
  1552.  
  1553. GEN padicff(f,p,r)
  1554.      GEN f,p;
  1555.      long r;
  1556.  
  1557. /* interne donc aucune verification */
  1558. /* En entree, res est un polynome a coeffs dans Z_p sans facteur carre */
  1559.  
  1560. {
  1561.   long av=avma,tetpil,lx,n,i,j,k,v,fl2;
  1562.   GEN rac,pro,p1,p2,y,yi,quatre;
  1563.  
  1564.   fl2=gegal(p,gdeux);
  1565.   if(fl2&&(r>=2)) err(impl,"padicfactor for p=2");
  1566.   rac=factmod(f,p);lx=lg(rac[1]);
  1567.   if(r==1)
  1568.     {
  1569.       tetpil=avma;y=cgetg(lx,18);p2=(GEN)rac[1];
  1570.       for(i=1;i<lx;i++) y[i]=(long)gcvtop(p2[i],p,1);
  1571.       return gerepile(av,tetpil,y);
  1572.     }
  1573.   if(fl2) quatre=stoi(4);
  1574.   n=lgef(f)-3;pro=cgetg(n+1,18);j=0;v=varn(f);p2=lift(rac[1]);
  1575.   for(i=1;i<lx;i++)
  1576.     {
  1577.       p1=(GEN)p2[i];
  1578.       if(lgef(p1)==4)
  1579.     {
  1580.       p1=gcmp0(p1[2])?(GEN)p1[2]:gsub(fl2?quatre:p,p1[2]);
  1581.       yi=cgetg(5,7);yi[2]=lclone(p);
  1582.       if(signe(p1))
  1583.         {
  1584.           if((!fl2)||mpodd(p1))
  1585.         {
  1586.           yi[3]=lpuigs(p,r);yi[4]=(long)p1;yi[1]=0x8000+(r<<16);
  1587.         }
  1588.           else
  1589.         {
  1590.           yi[3]=lpuigs(p,r);yi[4]=un;yi[1]=0x8001+(r<<16);
  1591.         }
  1592.         }
  1593.       else {yi[3]=un;yi[4]=zero;yi[1]=0x8000+r;}
  1594.       p1=apprgen(f,yi);
  1595.       for(k=1;k<lg(p1);k++) pro[++j]=lsub(polx[v],p1[k]);
  1596.     }
  1597.       else
  1598.     {
  1599.       yi=gmodulcp(gcvtop(polx[v],p,r),gcvtop(p2[i],p,r));
  1600.       p1=apprgen9(f,yi);
  1601.       for(k=1;k<lg(p1);k++) pro[++j]=(long)caract(p1[k],v);
  1602.     }
  1603.     }
  1604.   tetpil=avma;p1=cgetg(j+1,18);for(i=1;i<=j;i++) p1[i]=lcopy(pro[i]);
  1605.   return gerepile(av,tetpil,p1);
  1606. }  
  1607.  
  1608. GEN factorpadic(x,p,r)
  1609.   GEN x,p;
  1610.   long r;
  1611.  
  1612. {
  1613.   long av=avma,av2,vv,k,i,j,i1,f,nbfac;
  1614.   GEN res,p1,p2,y,d,fa[100],a,ap,t,v,w;
  1615.   
  1616.   if((typ(x)!=10)||(!signe(x))) err(rootper1);
  1617.   if(r<=0) err(rootper4);
  1618.   y=cgetg(3,19);if(lgef(x)==3) {y[1]=lgetg(1,18);y[2]=lgetg(1,18);return y;}
  1619.   if(lgef(x)==4)
  1620.     {
  1621.       p1=cgetg(2,18);y[1]=(long)p1;p1[1]=lcopy(x);
  1622.       p1=cgetg(2,18);y[2]=(long)p1;p1[1]=un;return y;
  1623.     }
  1624.   d=content(x);vv=varn(x);a=gdiv(x,d);ap=deriv(a,vv);t=ggcd(a,ap);v=gdiv(a,t);
  1625.   w=gdiv(ap,t);j=0;f=1;nbfac=0;
  1626.   while(f)
  1627.     {
  1628.       j++;w=gsub(w,deriv(v,vv));f=signe(w);
  1629.       if(f)
  1630.     {
  1631.       res=ggcd(v,w);v=gdiv(v,res);w=gdiv(w,res);
  1632.     }
  1633.       else res=v;
  1634.       fa[j]=(lgef(res)>3) ? padicff(res,p,r) : cgetg(1,18);
  1635.       nbfac+=(lg(fa[j])-1);
  1636.     }
  1637.   av2=avma;y=cgetg(3,19);p1=cgetg(nbfac+1,18);y[1]=(long)p1;
  1638.   p2=cgetg(nbfac+1,18);y[2]=(long)p2;
  1639.   for(i=1,k=0;i<=j;i++)
  1640.     for(i1=1;i1<lg(fa[i]);i1++)
  1641.       {
  1642.     p1[++k]=lcopy(fa[i][i1]);p2[k]=lstoi(i);
  1643.       }
  1644.   return gerepile(av,av2,y);
  1645. }
  1646.  
  1647.  
  1648. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1649. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1650. /*                                                                 */
  1651. /*                     FACTORISATION DANS F_q                      */
  1652. /*                                                                 */
  1653. /*                                                                 */
  1654. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1655. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1656.  
  1657. GEN     factmod9(f,p,a)
  1658.      
  1659.      GEN     f,p,a;
  1660.      
  1661. {
  1662.   long  i,j,k,d,e,vf,va,expos[100],psim,nbfact,nbf,av,pk,tetpil,calc,qqs;
  1663.   GEN   y,t[100],f1,f2,f3,df1,df2,g,g1,g2;
  1664.   GEN   xmod,u,v,pd,q,qq,unfp,unfq;
  1665.   
  1666.   if((typ(a)!=10)||(typ(f)!=10)||gcmp0(f)||gcmp0(a)) err(factmoder);
  1667.   if(lgef(f)==3)
  1668.     {y=cgetg(3,19);y[1]=lgetg(1,18);y[2]=lgetg(1,18);return y;}
  1669.   vf=varn(f);va=varn(a);av=avma;unfp=gmodulcp(un,p);
  1670.   a=gmul(unfp,a);unfq=gmodulcp(gmul(unfp,polun[va]),a);
  1671.   f=gmul(unfq,f);qq=gpuigs(p,lgef(a)-3);qqs=itos(qq);
  1672.   
  1673.   /*  1ere etape : trouver les facteurs square-free produit
  1674.       des facteurs premiers de meme multiplicite    */
  1675.  
  1676.   f1=f;e=0;nbfact=1;pk=1;df1=deriv(f1,vf);calc=0;
  1677.  
  1678.   do
  1679.     {
  1680. /*printf("debut 1re etape\n");*/
  1681.       e+=pk;
  1682.       while (gcmp0(df1))
  1683.         {
  1684.           pk *= (psim=itos(p));e=pk;
  1685.           j=(lgef(f1)-3)/psim+3;f2=cgetg(j,10);
  1686.           f2[1]=0x1000000+j;setvarn(f2,vf);
  1687.           for(i=2;i<j;i++)
  1688.             f2[i]=lcopy(f1[psim*(i-2)+2]);
  1689.           f1=f2;df1=deriv(f1,vf);calc=0;
  1690.     }
  1691.       if(calc) f2=f3;else f2=ggcd(f1,df1);calc=1;
  1692.       if (lgef(f2)<4) u=f1;
  1693.       else
  1694.         {
  1695.           g1=gdiv(f1,f2);
  1696.           if (gcmp0(df2=deriv(f2,vf)))
  1697.             {
  1698.               u=g1;f3=f2;
  1699.         }
  1700.           else
  1701.             {
  1702.               f3=ggcd(f2,df2);
  1703.               if (lgef(f3)<4) u=gdiv(g1,f2);
  1704.               else
  1705.                 {
  1706.                   g2=gdiv(f2,f3);u=gdiv(g1,g2);
  1707.          }
  1708.         }
  1709.     }
  1710.       
  1711.       /*  2eme etape : 
  1712.           Ici u est un polynome square-free
  1713.           (produit des facteurs premiers de meme multiplicite  e :
  1714.           trouver les facteurs (square-free) produit des facteurs de meme
  1715.       degre d */
  1716.       
  1717.       d=0;pd=gun;
  1718.       xmod=gmodulcp(polx[vf],u);v=xmod;
  1719.       while(d<(lgef(u)-3)>>1)
  1720.         {
  1721.                                        /*printf("debut 2me etape\n");*/
  1722.           d++;
  1723.           pd=mulii(pd,qq);
  1724.           q=shifti(subis(pd,1),-1);
  1725.           v=gpui(v,qq);
  1726.           g=ggcd((gsub(v,xmod))[2],u);
  1727.       
  1728.           if (lgef(g)>3)
  1729.             {
  1730.                                       /*printf("debut 3me etape\n");*/
  1731.           
  1732. /*  3eme etape :
  1733. Ici g est produit de pol irreductibles ayant tous le meme degre d;*/
  1734.  
  1735.           t[nbfact]=g;j=nbfact+(lgef(g)-3)/d;
  1736.           split9(qqs,t+nbfact,d,p[2],q,unfq,qqs,a);
  1737. /* le premier parametre est un entier variable m qui sera converti en un
  1738.    polynome w dont les coeff sont ses digits en base p (initialement m = p --> X)
  1739.    pour faire pgcd de g avec w^(qq^d-1)/2 jusqu'a casser. */
  1740.           for(;nbfact<j;expos[nbfact++]=e);
  1741.           u=gdiv(u,g);v=gmodulcp(v[2],u);
  1742.         }
  1743.                                           /*printf("fin 3me etape\n");*/
  1744.  
  1745.     }                                 /*printf("fin 2me etape\n");*/
  1746.       if (lgef(u)>3) {t[nbfact]=u;expos[nbfact++]=e;}
  1747.       f1=f2;df1=df2;
  1748.     }
  1749.                                           /*printf("fin 1re etape\n");*/
  1750.   while(lgef(f1)>3);
  1751.   
  1752.                                     /*printf("fin de l'algorithme\n");*/
  1753.   nbf=nbfact;
  1754.   tetpil=avma;
  1755.   t[1]=gdiv(t[1],((GEN)t[1])[lgef(t[1])-1]);
  1756.   for(j=2;j<nbfact;j++)
  1757.     {
  1758.       if (expos[j]) t[j]=gdiv(t[j],((GEN)t[j])[lgef(t[j])-1]);
  1759.       for (k=1;k<j;k++)
  1760.     if (expos[k]&&polegal(t[j],t[k])) {expos[k]+=expos[j];expos[j]=0;nbf--;k=j;}
  1761.     }
  1762.   y=cgetg(3,19);
  1763.   u=cgetg(nbf,18);y[1]=(long)u;
  1764.   v=cgetg(nbf,18);y[2]=(long)v;
  1765.   for(j=1,k=0;j<nbfact;j++)
  1766.     {
  1767.       if (expos[j])
  1768.         {
  1769.           k++;
  1770.           u[k]=(long)t[j];
  1771.           v[k]=lstoi(expos[j]);
  1772.     }
  1773.     }
  1774.   return(gerepile(av,tetpil,y));
  1775. }
  1776.  
  1777. GEN stopoly9(m,p,qqs,v,a)
  1778.      long p,v,m,qqs;
  1779.      GEN a;
  1780.  
  1781. /* renvoie un polynome de la variable v
  1782. dont les coef sont les digits de m en base qq */
  1783.  
  1784. {
  1785.   GEN y,p2;
  1786.   long p1,l=2,l1,i,j,va=varn(a),c[1000],d[1000];
  1787.  
  1788.   do {c[l++]=m%qqs;m/=qqs;} while(m);
  1789.   y=cgetg(l,10);y[1]=0x1000000+l;setvarn(y,v);
  1790.   for(i=2;i<l;i++) 
  1791.     {
  1792.       p1=c[i];l1=2;
  1793.       do {d[l1++]=p1%p;p1/=p;} while(p1);
  1794.       p2=cgetg(l1,10);p2[1]=0x1000000+l1;setvarn(p2,va);
  1795.       for(j=2;j<l1;j++) p2[j]=lstoi(d[j]);
  1796.       y[i]=lmodulcp(p2,a);
  1797.     }
  1798.   return y;
  1799. }
  1800.  
  1801. GEN stopoly92(d1,v,a,ptres)
  1802.      long v,d1;
  1803.      GEN a,*ptres;
  1804.  
  1805. /* renvoie un polynome aleatoire de la variable v
  1806. de degre inferieur ou egal a 2*d1-1 */
  1807.  
  1808. {
  1809.   GEN y,p2;
  1810.   long p1,l1,i,j,d2,l,va=varn(a),c[1000],d[1000],k=lgef(a)-3;
  1811.  
  1812. /* qqs=2^k */
  1813.  
  1814.   c[1]=1;d2=d1+d1+1;
  1815.   do
  1816.     {
  1817.       for(l=2;l<=d2;l++) c[l]=(rand())>>(31-k);
  1818.       l=d2;while(!c[l]) l--;
  1819.     }
  1820.   while(l<=2);
  1821.   l++;
  1822.   y=cgetg(l,10);y[1]=0x1000000+l;setvarn(y,v);
  1823.   for(i=2;i<l;i++) 
  1824.     {
  1825.       p1=c[i];l1=2;
  1826.       do {d[l1++]=p1&1;p1>>=1;} while(p1);
  1827.       p2=cgetg(l1,10);p2[1]=0x1000000+l1;setvarn(p2,va);
  1828.       for(j=2;j<l1;j++) p2[j]=lstoi(d[j]);
  1829.       y[i]=lmodulcp(p2,a);
  1830.     }
  1831.   p1=(rand())>>(31-k);l1=2;
  1832.   do {d[l1++]=p1&1;p1>>=1;} while(p1);
  1833.   p2=cgetg(l1,10);p2[1]=0x1000000+l1;setvarn(p2,va);
  1834.   for(j=2;j<l1;j++) p2[j]=lstoi(d[j]);
  1835.   *ptres=gmodulcp(p2,a);
  1836.   return y;
  1837. }
  1838.  
  1839. split9(m,t,d,p,q,unfq,qqs,a)
  1840.      GEN *t,q,unfq,a;
  1841.      long m,d,p,qqs;
  1842.  
  1843. /*----------------------------------------------------- 
  1844.    Programme recursif :
  1845.   Entree:
  1846.    m entier arbitraire ( converti en un polynome w )
  1847.    p nb premier; q=(qq^d-1)/2
  1848.    t[0] polynome de degre k*d prod de k fact de deg d.
  1849.   Sortie:
  1850.    t[0],t[1]...t[k-1] contiennent les k facteurs de g
  1851. ------------------------------------------------------*/
  1852.  
  1853. {  
  1854.   long j,l,v,av,av1,dv;
  1855.   GEN w,w0,wm,res;
  1856.  
  1857.   if ((dv=lgef(*t)-3)==d) return 0;
  1858.   v=varn(*t);
  1859.   do
  1860.     {
  1861.       av=avma;
  1862.       if(p==2)
  1863.     {
  1864.       w=gmul(stopoly92(d,v,a,&res),unfq);
  1865.       for(w0=w,j=1;j<d;j++) w=gmod(gadd(w0,gpuigs(w,qqs)),*t);
  1866.       w=gsub(w,res);
  1867.     }
  1868.       else
  1869.     {
  1870.       w=gmul(stopoly9(m,p,qqs,v,a),unfq);m++;
  1871.       wm=gpui(gmodulcp(w,*t),q);w=gsub(wm[2],unfq);
  1872.     }
  1873.       av1=avma;
  1874.       w=gerepile(av,av1,ggcd(*t,w));
  1875.       l=lgef(w)-3;
  1876.     }
  1877.   while((!l)||(l==dv));
  1878.   t[j=l/d]=gdiv(*t,w);*t=w;
  1879.   split9(m,t+j,d,p,q,unfq,qqs,a);split9(m,t,d,p,q,unfq,qqs,a);
  1880. }
  1881.  
  1882. GEN randompoly9(d1,p,v,a)
  1883.      long v,d1;
  1884.      GEN a,p;
  1885.  
  1886. /* renvoie un polynome aleatoire de la variable v
  1887. de degre inferieur ou egal a d1 */
  1888.  
  1889. {
  1890.   GEN y,p2;
  1891.   long p1,ps,qqs,l1,i,j,l,va=varn(a),c[1000],d[1000],k=lgef(a)-3;
  1892.  
  1893.   c[1]=1;qqs=itos(gpuigs(p,k));ps=itos(p);
  1894.   do
  1895.     {
  1896.       for(l=2;l<=d1+2;l++) c[l]=(rand()*(((double)qqs)/C31));
  1897.       l=d1+2;while(!c[l]) l--;
  1898.     }
  1899.   while(l<=2);
  1900.   l++;
  1901.   y=cgetg(l,10);y[1]=0x1000000+l;setvarn(y,v);
  1902.   for(i=2;i<l;i++) 
  1903.     {
  1904.       p1=c[i];l1=2;
  1905.       do {d[l1++]=p1%ps;p1/=ps;} while(p1);
  1906.       p2=cgetg(l1,10);p2[1]=0x1000000+l1;setvarn(p2,va);
  1907.       for(j=2;j<l1;j++) p2[j]=lstoi(d[j]);
  1908.       y[i]=lmodulcp(p2,a);
  1909.     }
  1910.   return y;
  1911. }
  1912.  
  1913.  
  1914.