home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Professional / OS2PRO194.ISO / os2 / progs / pari / pari_137 / src / gen3.c < prev    next >
Text File  |  1992-05-20  |  52KB  |  2,148 lines
  1. /********************************************************************/
  2. /********************************************************************/
  3. /**                                                                **/
  4. /**                +++++++++++++++++++++++++++++++                 **/
  5. /**                +                             +                 **/
  6. /**                +    OPERATIONS GENERIQUES    +                 **/
  7. /**                +     (troisieme partie)      +                 **/
  8. /**                +                             +                 **/
  9. /**                +     copyright Babe Cool     +                 **/
  10. /**                +                             +                 **/
  11. /**                +++++++++++++++++++++++++++++++                 **/
  12. /**                                                                **/
  13. /**                                                                **/
  14. /********************************************************************/
  15. /********************************************************************/
  16.  
  17. # include "genpari.h"
  18.  
  19. /*******************************************************************/
  20. /*******************************************************************/
  21. /*                                                                 */
  22. /*                 LISTE DES TYPES GENERIQUES                      */
  23. /*                 ~~~~~~~~~~~~~~~~~~~~~~~~~~                      */
  24. /*                                                                 */
  25. /*  1  :entier long     [ cod1 ] [ cod2 ] [ man1 ] ... [ manl ]    */
  26. /*  2  :reel            [ cod1 ] [ cod2 ] [ man1 ] ... [ manl ]    */
  27. /*  3  :entier modulo   [ code ] [ mod  ] [ entier modulo ]        */
  28. /*  4  :fraction        [ code ] [ num. ] [ den. ]                 */
  29. /*  5  :nfraction       [ code ] [ num. ] [ den. ]                 */
  30. /*  6  :complexe        [ code ] [ reel ] [ imag ]                 */
  31. /*  7  :p-adique        [ cod1 ] [ cod2 ] [ p ] [ p^r ] [ entier]  */
  32. /*  8  :quadrat         [ cod1 ] [ mod  ] [ reel ] [ imag ]        */
  33. /*  9  :poly mod        [ code ] [ mod  ] [ polynome  mod ]        */
  34. /* --------------------------------------------------------------- */
  35. /*  10 :polynome        [ cod1 ] [ cod2 ] [ man1 ] ... [ manl ]    */
  36. /*  11 :serie           [ cod1 ] [ cod2 ] [ man1 ] ... [ manl ]    */
  37. /*  13 :fr.rat          [ code ] [ num. ] [ den. ]                 */
  38. /*  14 :n.fr.rat        [ code ] [ num. ] [ den. ]                 */
  39. /*  15 :formqre         [ code ] [  a  ] [  b  ] [  c  ] [ del ]   */
  40. /*  16 :formqim         [ code ] [  a   ] [  b   ] [  c   ]        */
  41. /*  17 :vecteur ligne   [ code ] [  x1  ] ... [  xl  ]             */
  42. /*  18 :vecteur colonne [ code ] [  x1  ] ... [  xl  ]             */
  43. /*  19 :matrice         [ code ] [ col1 ] ... [ coll ]             */
  44. /*                                                                 */
  45. /*******************************************************************/
  46. /*******************************************************************/
  47.  
  48. /********************************************************************/
  49. /********************************************************************/
  50. /**                                                                **/
  51. /**                        TYPE D'UN GEN                           **/
  52. /**                                                                **/
  53. /**                    (POUR GP UNIQUEMENT)                        **/
  54. /**                                                                **/
  55. /********************************************************************/
  56. /********************************************************************/
  57.  
  58. GEN gtype(x)
  59.   
  60.    GEN  x;
  61.   
  62. {
  63.   return stoi(typ(x));
  64. }
  65.  
  66. /********************************************************************/
  67. /********************************************************************/
  68. /**                                                                **/
  69. /**                 NUMERO DE LA VARIABLE PRINCIPALE               **/
  70. /**                                                                **/
  71. /********************************************************************/
  72. /********************************************************************/
  73.  
  74. int    gvar(x)
  75.   
  76.      GEN  x;
  77.   
  78. {
  79.   long tx=typ(x),i,v,w;
  80.   if((tx>>1)==5) return varn(x);
  81.   if((tx<9)||(tx==15)||(tx==16)) return BIGINT;
  82.   if(tx==9) return varn(x[1]);
  83.   v=BIGINT;
  84.   for(i=1;i<lg(x);i++) {w=gvar(x[i]);if(w<v) v=w;}
  85.   return v;
  86. }
  87.  
  88. int gvar2(x)
  89.      GEN x;
  90.  
  91. {
  92.   long tx=typ(x),i,v,w;
  93.   if((tx<9)||(tx==15)||(tx==16)) return BIGINT;
  94.   if(tx==9) {v=gvar2(x[1]);w=gvar2(x[2]);return min(v,w);}
  95.   v=BIGINT;
  96.   if((tx==11)&&(!signe(x))) return v;
  97.   if(tx<12) 
  98.     {
  99.       for(i=2;i<((tx==11)?lg(x):lgef(x));i++)
  100.     {w=gvar(x[i]);if(w<v) v=w;} 
  101.       return v;
  102.     }
  103.   else {for(i=1;i<lg(x);i++) {w=gvar2(x[i]);if(w<v) v=w;} return v;}
  104. }
  105.  
  106. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  107. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  108. /*                                                                 */
  109. /*            PRECISION D'UN SCALAIRE                              */
  110. /*                                                                 */
  111. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  112. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  113.  
  114. int precision(x)
  115.      
  116.      GEN     x;
  117.      
  118. {
  119.   long    l1,l2,tx=typ(x);
  120.   
  121.   if (tx==2) return max(lg(x),2-(expo(x)>>5));
  122.   if (tx==6)
  123.   {
  124.     l1=precision(x[1]);l2=precision(x[2]);
  125.     if(!l1) return l2;
  126.     if(!l2) return l1;
  127.     return (l1>l2) ? l2 : l1;
  128.   }
  129.   return 0;
  130. }
  131.  
  132. /* degre de x par rapport a la variable principale */
  133. /* et 0 si x est nul */
  134.  
  135. int tdeg(x)
  136.      GEN x;
  137.      
  138. {
  139.   long tx=typ(x);
  140.   
  141.   if(gcmp0(x)) return 0;
  142.   if(tx<10) return 0;
  143.   switch(tx)
  144.   {
  145.     case 10: return lgef(x)-3;
  146.     case 13:
  147.     case 14: return tdeg(x[1])-tdeg(x[2]);
  148.     default: err(tdeger);
  149.   }
  150. }
  151.  
  152. /********************************************************************/
  153. /********************************************************************/
  154. /**                                                                **/
  155. /**                     MULTIPLICATION SIMPLE                      **/
  156. /**                                                                **/
  157. /********************************************************************/
  158. /********************************************************************/
  159.  
  160. GEN   gmulsg(s,y)
  161.   
  162.    long  s;
  163.    GEN   y;
  164.   
  165. {
  166.   long ty=typ(y),ly=lg(y),i,l,tetpil;
  167.   GEN  z,p1;
  168.  
  169.   switch(ty)
  170.   {
  171.   case 1 : z=mulsi(s,y);break;
  172.   
  173.   case 2 : z=mulsr(s,y);break;
  174.   
  175.   case 3 : z=cgetg(ly,ty);z[1]=copyifstack(y[1]);
  176.     l=avma;
  177.     p1=mulsi(s,y[2]);
  178.     tetpil=avma;
  179.     z[2]=lpile(l,tetpil,modii(p1,y[1]));
  180.     break;
  181.   
  182.   case 4 :
  183.   
  184.   case 5 : z=cgetg(ly,ty);
  185.     z[1]=lmulsi(s,y[1]);
  186.     z[2]=lcopy(y[2]);
  187.     if (ty==4) gredsp(&z);
  188.     break;
  189.   
  190.   case 6 : z=cgetg(ly,ty);
  191.     z[1]=lmulsg(s,y[1]);
  192.     z[2]=lmulsg(s,y[2]);
  193.     break;
  194.   
  195.   case 7 : if(s)
  196.     {
  197.       l=avma;p1=cgetp(y);gaffsg(s,p1);
  198.       tetpil=avma;z=gerepile(l,tetpil,gmul(p1,y));
  199.     }
  200.     else z=gzero;
  201.     break;
  202.   
  203.   case 8 : z=cgetg(ly,ty);
  204.     z[2]=lmulsg(s,y[2]);
  205.     z[3]=lmulsg(s,y[3]);
  206.     z[1]=copyifstack(y[1]);
  207.     break;
  208.  
  209.   case 9 : z=cgetg(ly,ty);z[2]=lmulsg(s,y[2]);
  210.     z[1]=copyifstack(y[1]);
  211.     break;
  212.   
  213.   case 10: 
  214.     if ((!s)||gcmp0(y)) {z=cgetg(2,10);z[1]=2;setvarn(z,varn(y));}
  215.     else
  216.     {
  217.       ly=lgef(y);z=cgetg(ly,ty);
  218.       for (i=2;i<ly;i++)
  219.       z[i]=lmulsg(s,y[i]);
  220.       z[1]=y[1];normalizepol(&z);
  221.     }
  222.     break;
  223.   
  224.   case 11: if (!s)
  225.     {
  226.     z=cgetg(3,10);z[1]=2;setvarn(z,varn(y));
  227.     }
  228.   else
  229.     {
  230.     if (gcmp0(y)) z=gcopy(y);
  231.     else
  232.       {
  233.       z=cgetg(ly,ty);
  234.       for (i=2;i<ly;i++)
  235.         z[i]=lmulsg(s,y[i]);
  236.       z[1]=y[1];
  237.       normalize(&z);
  238.       }
  239.     }
  240.     break;
  241.   
  242.   case 13:
  243.     if(s) 
  244.       {
  245.     l=avma;z=cgetg(ly,ty);z[1]=lmulsg(s,y[1]);z[2]=y[2];
  246.     tetpil=avma;z=gerepile(l,tetpil,gred(z));
  247.       }
  248.     else {z=cgetg(2,10);z[1]=2;setvarn(z,gvar(y));}
  249.     break;
  250.   case 14: 
  251.     if(s) {z=cgetg(ly,ty);z[1]=lmulsg(s,y[1]);z[2]=lcopy(y[2]);}
  252.     else {z=cgetg(2,10);z[1]=2;setvarn(z,gvar(y));}
  253.     break;
  254.   
  255.   case 17:
  256.   case 18:
  257.   case 19: z=cgetg(ly,ty);
  258.     for (i=1;i<ly;i++)
  259.     z[i]=lmulsg(s,y[i]);
  260.     break;
  261.   
  262.   default: err(gmuler1);
  263.   
  264.   }
  265.   return z;
  266. }
  267.  
  268. /********************************************************************/
  269. /********************************************************************/
  270. /**                                                                **/
  271. /**                       DIVISION SIMPLE                          **/
  272. /**                                                                **/
  273. /********************************************************************/
  274. /********************************************************************/
  275.  
  276. GEN   gdivgs(x,s)
  277.   
  278.    GEN   x;
  279.    long  s;
  280.   
  281. {
  282.   long  tx=typ(x),lx=lg(x),i,l,tetpil,court[3];
  283.   GEN z,p1;
  284.  
  285.   if(!s) err(gdiver2);
  286.   switch(tx)
  287.   {
  288.     case 1 : l=avma;z=dvmdis(x,s,&p1);
  289.       if(!signe(p1)) cgiv(p1);
  290.       else
  291.       {
  292.         avma=l;
  293.         z=cgetg(3,4);z[1]=lcopy(x);
  294.         z[2]=lstoi(s);
  295.         if(s>=0)
  296.         {
  297.           z[1]=lcopy(x);z[2]=lstoi(s);
  298.         }
  299.         else
  300.         {
  301.           z[1]=lnegi(x);z[2]=lstoi(-s);
  302.         }
  303.         gredsp(&z);
  304.       }
  305.       break;
  306.     
  307.     case 2 : z=divrs(x,s);break;
  308.     
  309.     case 4 :
  310.     case 5 : z=cgetg(lx,tx);z[1]=lcopy(x[1]);
  311.       z[2]=lmulsg(s,x[2]);
  312.       if(signe(z[2])<0)
  313.       {
  314.         mpnegz(z[1],z[1]);
  315.         mpnegz(z[2],z[2]);
  316.       }
  317.       if (tx==4) gredsp(&z);
  318.       break;
  319.     
  320.     case 6 : z=cgetg(lx,tx);
  321.       z[1]=ldivgs(x[1],s);
  322.       z[2]=ldivgs(x[2],s);
  323.       break;
  324.     
  325.     case 8 : z=cgetg(lx,tx);
  326.       z[1]=copyifstack(x[1]);
  327.       for (i=2;i<4;i++)
  328.       z[i]=ldivgs(x[i],s);
  329.       break;
  330.     
  331.     case 9 : z=cgetg(lx,tx);z[2]=ldivgs(x[2],s);
  332.       z[1]=copyifstack(x[1]);
  333.       break;
  334.     
  335.     case 10: lx=lgef(x);z=cgetg(lx,tx);
  336.       for (i=2;i<lx;i++)
  337.       z[i]=ldivgs(x[i],s);
  338.       z[1]=x[1];
  339.       break;
  340.     
  341.     case 11:
  342.       if (gcmp0(x)) z=gcopy(x);
  343.       else
  344.       {
  345.         z=cgetg(lx,tx);
  346.         for (i=2;i<lx;i++)
  347.         z[i]=ldivgs(x[i],s);
  348.         z[1]=x[1];
  349.         normalize(&z);
  350.       }
  351.       break;
  352.     
  353.     case 13: l=avma;z=cgetg(lx,tx);z[1]=x[1];z[2]=lmulsg(s,x[2]);
  354.       tetpil=avma;z=gerepile(l,tetpil,gred(z));
  355.       break;
  356.     case 14: z=cgetg(lx,tx);z[1]=lcopy(x[1]);z[2]=lmulsg(s,x[2]);
  357.       break;
  358.     
  359.     case 17:
  360.     case 18:
  361.     case 19: z=cgetg(lx,tx);
  362.       for (i=1;i<lx;i++)
  363.       z[i]=ldivgs(x[i],s);
  364.       break;
  365.     
  366.     default: court[0] = 0x1010003; affsi(s,court);z=gdiv(x,court);
  367.   }
  368.   return z;
  369. }
  370.  
  371. GEN   gaddpex(x,y)
  372.   
  373.    /* addition d'un type entier ou rationnel avec un p-adique */
  374.    /* x doit etre entier ou rationnel et y   p-adique     */
  375.    /* a usage interne donc aucune verification de type.   */
  376.   
  377.    GEN   x,y;
  378.   
  379. {
  380.   GEN   z,p,p1,p2,p3;
  381.   long  e1,e2,e3,av,tetpil;
  382.  
  383.   if(gcmp0(x)) z=gcopy(y);
  384.   else
  385.   {
  386.     av=avma;z=cgetg(5,7);e1=valp(y);
  387.     p=(GEN)(y[2]);
  388.     if(typ(x)>=4)
  389.     {
  390.       e3=pvaluation(x[1],p,&p2);
  391.       e3-=pvaluation(x[2],p,&p3);
  392.       p2=gdiv(p2,p3);
  393.     }
  394.     else e3=pvaluation(x,p,&p2);
  395.     e2=signe(y[4])?e1+precp(y)-e3:e1-e3;
  396.     z[2]=(long)p;
  397.     if(e2<=0) {z[3]=un;setprecp(z,0);}
  398.     else
  399.     {
  400.       setprecp(z,e2);
  401.       if(e1-e3)
  402.       {
  403.         p1=gpuigs(p,e1-e3);
  404.         z[3]=lmul(y[3],p1);
  405.       }
  406.       else z[3]=lcopy(y[3]);
  407.     }
  408.     z[4]=lmod(p2,z[3]);setvalp(z,e3);
  409.     tetpil=avma;z=gerepile(av,tetpil,gadd(z,y));
  410.   }
  411.   return z;
  412. }
  413.  
  414. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  415. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  416. /*                                                                 */
  417. /*                    MODULO GENERAL                               */
  418. /*                                                                 */
  419. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  420. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  421.  
  422.  
  423. GEN     gmod(x,y)
  424.      
  425.      GEN     x,y;
  426.      
  427. {
  428.   long  tx,ty,l,tetpil,lx,i;
  429.   GEN   z,p1;
  430.   
  431.   ty=typ(y);tx=typ(x);
  432.   if(ty==1)
  433.     switch(tx)
  434.     {
  435.       case 1 : z=modii(x,y);break;
  436.       case 2 : err(gmoder1);break;
  437.       case 3 : z=cgetg(3,3);z[1]=lmppgcd(x[1],y);
  438.         z[2]=lmodii(x[2],z[1]);
  439.         break;
  440.       case 4 :
  441.       case 5 : l=avma;if(tx==5) p1=gred(x);else p1=x;
  442.         p1=mulii(p1[1],mpinvmod(p1[2],y));
  443.         tetpil=avma;
  444.         z=gerepile(l,tetpil,modii(p1,y));
  445.         break;
  446.       case 8 : z=cgetg(4,8);z[1]=copyifstack(x[1]);
  447.         z[2]=lmod(x[2],y);
  448.         z[3]=lmod(x[3],y);
  449.         break;
  450.       case 10: z=gzero;break;
  451.       case 17:
  452.       case 18:
  453.       case 19: lx=lg(x);z=cgetg(lx,tx);
  454.         for(i=1;i<lx;i++) z[i]=lmod(x[i],y);
  455.         break;
  456.       default: err(gmoder1);
  457.     }
  458.   else
  459.   {
  460.     if(ty==10)
  461.       if(tx<9) z=gcopy(x);
  462.       else
  463.       {
  464.         switch(tx)
  465.         {
  466.           case 9 : z=cgetg(3,9);z[1]=lgcd(x[1],y);
  467.             z[2]=lres(x[2],z[1]);
  468.             break;
  469.           case 10: z=gres(x,y);break;
  470.           case 11: err(gmoder3);break;
  471.           case 13:
  472.           case 14: l=avma;if(tx==14) p1=gred(x);else p1=x;
  473.             p1=gmul(p1[1],ginvmod(p1[2],y));
  474.             tetpil=avma;z=gerepile(l,tetpil,gres(p1,y));
  475.             break;
  476.           case 17:
  477.           case 18:
  478.           case 19: lx=lg(x);z=cgetg(lx,tx);
  479.             for(i=1;i<lx;i++) z[i]=lmod(x[i],y);
  480.             break;
  481.           default: err(gmoder3);
  482.         }
  483.       }
  484.     else err(gmoder5);
  485.   }
  486.   return z;
  487. }
  488.  
  489. GEN     gmodulo(x,y)
  490.      
  491.      GEN        x,y;
  492.      
  493. {
  494.   long    ty=typ(y);
  495.   GEN   z;
  496.   
  497.   if(ty==1)
  498.   {
  499.     if(typ(x)!=1) err(gmoder1);
  500.     z=cgetg(3,3);z[1] = lclone(y);
  501.     z[2]=lmodii(x,y);
  502.   }
  503.   else
  504.     if(ty==10)
  505.     {
  506.       z=cgetg(3,9);z[1] = lclone(y);
  507.       ty=typ(x);
  508.       if(ty>=10)
  509.     {if(ty==10) z[2]=lmod(x,y);else err(gmoder1);}
  510.       else z[2]=lcopy(x);
  511.     }
  512.     else err(gmoder1);
  513.   return z;
  514. }
  515.  
  516. GEN     gmodulcp(x,y)
  517.      
  518.      GEN        x,y;
  519.      
  520. {
  521.   long    ty=typ(y);
  522.   GEN   z;
  523.   
  524.   if(ty==1)
  525.   {
  526.     if(typ(x)!=1) err(gmoder1);
  527.     z=cgetg(3,3);z[1]=lcopy(y);
  528.     z[2]=lmodii(x,y);
  529.   }
  530.   else
  531.     if(ty==10)
  532.     {
  533.       z=cgetg(3,9);z[1]=lcopy(y);
  534.       ty=typ(x);
  535.       if(ty>=10)
  536.     {if(ty==10) z[2]=lmod(x,y);else err(gmoder1);}
  537.       else z[2]=lcopy(x);
  538.     }
  539.     else err(gmoder1);
  540.   return z;
  541. }
  542.  
  543.  
  544. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  545. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  546. /*                                                                 */
  547. /*                 DIVISION ENTIERE GENERALE                       */
  548. /*            DIVISION ENTIERE AVEC RESTE GENERALE                 */
  549. /*                                                                 */
  550. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  551. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  552.  
  553. GEN     gdivent(x,y)
  554.      
  555.      GEN     x,y;
  556.      
  557. {
  558.   long    tx=typ(x),ty=typ(y),av,tetpil,i;
  559.   GEN     z,p1;
  560.   
  561.   if (tx==1)
  562.   {
  563.     if(ty==1)
  564.     {
  565.       av=avma;z=dvmdii(x,y,&p1);
  566.       i=signe(p1);cgiv(p1);
  567.       if(i<0)
  568.       {
  569.         tetpil=avma;z=gerepile(av,tetpil,gaddgs(z,-signe(y)));
  570.       }
  571.     }
  572.     else
  573.     {
  574.       if(ty!=10) err(gdiventer);
  575.       z=gzero;
  576.     }
  577.     return z;
  578.   }
  579.   if((ty!=10)||(tx>10)) err(gdiventer);
  580.   if(tx==10) return gdeuc(x,y);
  581.   else return gzero;
  582. }
  583.  
  584. GEN     gdiventres(x,y)
  585.      
  586.      GEN     x,y;
  587.      
  588. {
  589.   long    tx=typ(x),ty=typ(y);
  590.   GEN     z;
  591.   
  592.   z=cgetg(3,18);
  593.   if (tx==1)
  594.   {
  595.     if(ty==1) z[1]=(long)dvmdii(x,y,z+2);
  596.     else
  597.     {
  598.       if(ty!=10) err(gdiventer);
  599.       z[1]=zero;z[2]=lcopy(x);
  600.     }
  601.   }
  602.   else
  603.   {
  604.     if((ty!=10)||(tx>10)) err(gdiventer);
  605.     if(tx==10) z[1]=ldivres(x,y,z+2);
  606.     else {z[1]=zero;z[2]=lcopy(y);}
  607.   }
  608.   return z;
  609. }
  610.  
  611. GEN     gdivmod(x,y,pr)
  612.      
  613.      GEN     x,y,*pr;
  614.      
  615. {
  616.   long    tx=typ(x),ty=typ(y);
  617.   
  618.   if(tx==1)
  619.   {
  620.     if(ty==1) return dvmdii(x,y,pr);
  621.     if(ty==10) {*pr=gcopy(x);return gzero;}
  622.     else err(gdivmoder);
  623.   }
  624.   if (tx==10) return poldivres(x,y,pr);
  625.   else err(gdivmoder);
  626. }
  627.  
  628.  
  629. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  630. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  631. /*                                                                 */
  632. /*                       SHIFT D'UN GEN                            */
  633. /*                                                                 */
  634. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  635. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  636.  
  637. /*      SHIFT TRONQUE SI n<0 (MULTIPLICATION TRONQUEE PAR 2**n)  */
  638.  
  639. GEN     gshift(x,n)
  640.      
  641.      GEN     x;
  642.      long    n;
  643.      
  644. {
  645.   long  tx,lx,i,l;
  646.   GEN   y;
  647.   lx=lg(x);tx=typ(x);
  648.   if(gcmp0(x)) y=gcopy(x);
  649.   else
  650.     switch(tx)
  651.     {
  652.       case 1 :
  653.       case 2 : y=mpshift(x,n);break;
  654.       case 17:
  655.       case 18:
  656.       case 19: y=cgetg(lx,tx);l=lontyp[tx];
  657.         for(i=1;i<l;y[i]=x[i],i++);
  658.         for(i=l;i<lx;i++)
  659.           y[i]=lshift(x[i],n);
  660.         break;
  661.       default: y=gmul2n(x,n);
  662.     }
  663.   return y;
  664. }
  665.  
  666. /*      SHIFT VRAI (MULTIPLICATION EXACTE PAR 2**n)     */
  667.  
  668. GEN     gmul2n(x,n)
  669.      
  670.      GEN     x;
  671.      long    n;
  672.      
  673. {
  674.   long  tx,lx,i,l,tetpil;
  675.   GEN   y,p1;
  676.   lx=lg(x);tx=typ(x);
  677.   if(gcmp0(x)) y=gcopy(x);
  678.   else
  679.     switch(tx)
  680.     {
  681.       case 1 : if(n>=0) y=mpshift(x,n);
  682.         else
  683.         {
  684.           y=cgetg(3,4);y[1]=lcopy(x);
  685.           y[2]=lmpshift(un,-n);gredsp(&y);
  686.         }
  687.         break;
  688.         
  689.       case 2 : y=mpshift(x,n);break;
  690.       case 4 :
  691.       case 5 : y=cgetg(lx,tx);
  692.         if(n>=0)
  693.         {
  694.           y[1]=lmpshift(x[1],n);
  695.           y[2]=lcopy(x[2]);
  696.         }
  697.         else
  698.         {
  699.           y[2]=lmpshift(x[2],-n);
  700.           y[1]=lcopy(x[1]);
  701.         }
  702.         if(tx==4) gredsp(&y);
  703.         break;
  704.       case 8 : y=cgetg(lx,tx);
  705.     y[1]=copyifstack(x[1]);
  706.         for(i=2;i<lx;i++)
  707.           y[i]=lmul2n(x[i],n);
  708.         break;
  709.       case 9 : y=cgetg(lx,tx);
  710.     y[1]=copyifstack(x[1]);
  711.     y[2]=lmul2n(x[2],n);
  712.     break;
  713.  
  714.       case 10: lx=lgef(x);
  715.       case 6 :
  716.       case 11:
  717.       case 17:
  718.       case 18:
  719.       case 19: y=cgetg(lx,tx);l=lontyp[tx];
  720.         for(i=1;i<l;y[i]=x[i],i++);
  721.         for(i=l;i<lx;i++)
  722.           y[i]=lmul2n(x[i],n);
  723.         break;
  724.       case 13: l=avma;y=cgetg(lx,tx);
  725.         if(n>=0) {y[1]=lmul2n(x[1],n);y[2]=x[2];}
  726.         else {y[2]=lmul2n(x[2],-n);y[1]=x[1];}
  727.     tetpil=avma;y=gerepile(l,tetpil,gred(y));
  728.         break;
  729.       case 14: y=cgetg(lx,tx);
  730.     if(n>=0) {y[1]=lmul2n(x[1],n);y[2]=lcopy(x[2]);}
  731.         else {y[2]=lmul2n(x[2],-n);y[1]=lcopy(x[1]);}
  732.         break;
  733.       case 3 :
  734.       case 7 : l=avma;p1=gmul2n(un,n);tetpil=avma;
  735.         y=gerepile(l,tetpil,gmul(p1,x));
  736.         break;
  737.       default: err(gmul2ner1);
  738.     }
  739.   return y;
  740. }
  741.  
  742.  
  743. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  744. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  745. /*                                                                 */
  746. /*                    INVERSE D' UN GEN                            */
  747. /*                                                                 */
  748. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  749. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  750.  
  751. GEN     ginv(x)
  752.      
  753.      GEN     x;
  754.      
  755. {
  756.   long    tx=typ(x),av,tetpil;
  757.   GEN     y;
  758.  
  759.   if(tx==16) 
  760.     {
  761.       y=gcopy(x);
  762.       if(cmpii(x[1],x[2])&&cmpii(x[1],x[3])) setsigne(y[2],-signe(y[2]));
  763.       return y;
  764.     }
  765.   if(tx==15) 
  766.     {
  767.       av=avma;y=gcopy(x);setsigne(y[2],-signe(y[2]));
  768.       setsigne(y[4],-signe(y[4]));tetpil=avma;
  769.       return gerepile(av,tetpil,redreal(y));
  770.     }
  771.   if (tx<15) return gdivsg(1,x);
  772.   if (tx==19) return invmat(x);
  773.   err(ginver);
  774. }
  775.  
  776.  
  777. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  778. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  779. /*                                                                 */
  780. /*           SUBSTITUTION DANS UN POLYNOME OU UNE SERIE            */
  781. /*                                                                 */
  782. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  783. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  784.  
  785. GEN     gsubst(x,v,y)
  786.      
  787.      GEN     x,y;
  788.      long    v;
  789.      
  790. {
  791.   long  tx,ty,l,lx,ly,vx,vy,e,ex,ey,tetpil;
  792.   long  av,av1,av2,av3,av4,av5,av6,av7,i,j,k,jb,decal;
  793.   GEN   t,p1,p2,z;
  794.   
  795.   tx=typ(x);ty=typ(y);lx=lg(x);ly=lg(y);
  796.   if ((ty>=15)&&(ty<=18)) err(gsubser2);
  797.   if ((ty==19)&&(ly!=lg(y[1]))) err(gsubser3);
  798.   if ((tx<9)||((tx==9)&&(v<=varn(x[1]))))
  799.     if(ty==19) z=gscalmat(x,ly-1);else z=gcopy(x);
  800.   else switch(tx)
  801.   {
  802.     case 10:
  803.       if(!signe(x))
  804.       {
  805.         if(ty==19) z=gscalmat(zero,ly-1);else z=gzero;
  806.       }
  807.       else
  808.       {
  809.     vx=varn(x);
  810.         if((ty==6)&&(v==vx)) z=poleval(x,y);
  811.         else
  812.         {
  813.           l=lgef(x);
  814.           if(ty!=19)
  815.           {
  816.             if(vx>v) z=gcopy(x);
  817.             else
  818.             {
  819.               if(vx==v)
  820.               {
  821.                 if(l==3) z=gcopy(x[2]);
  822.                 else
  823.                 {
  824.                   av=avma;z=(GEN)(x[l-1]);
  825.                   for (i=l-1;i>3;i--)
  826.                     z=gadd(gmul(z,y),x[i-1]);
  827.                   z=gmul(z,y);tetpil=avma;
  828.                   z=gerepile(av,tetpil,gadd(z,x[2]));
  829.                 }
  830.               }
  831.               else
  832.               {
  833.                 if(l==3) z=gsubst(x[2],v,y);
  834.                 else
  835.                 {
  836.                   av=avma;p1=polx[vx];
  837.                   z= gsubst(x[l-1],v,y);
  838.                   for (i=l-1;i>3;i--)
  839.                     z=gadd(gmul(z,p1),gsubst(x[i-1],v,y));
  840.                   z=gmul(z,p1);p2=gsubst(x[2],v,y);
  841.                   tetpil=avma;
  842.                   z=gerepile(av,tetpil,gadd(z,p2));
  843.                 }
  844.               }
  845.             }
  846.           }
  847.           else
  848.           {
  849.             if(ly!=lg(y[1])) err(gsubser3);
  850.             if(vx>v) z=gscalmat(x,ly-1);
  851.             else
  852.             {
  853.               if(vx==v)
  854.               {
  855.                 if(l==3) z=gscalmat(x[2],ly-1);
  856.                 else
  857.                 {
  858.                   av=avma;z=(GEN)(x[l-1]);
  859.                   for (i=l-1;i>3;i--)
  860.                     z=gaddmat(x[i-1],gmul(z,y));
  861.                   z=gmul(z,y);tetpil=avma;
  862.                   z=gerepile(av,tetpil,gaddmat(x[2],z));
  863.                 }
  864.               }
  865.               else
  866.               {
  867.                 if(l==3) z=gsubst(x[2],v,y);
  868.                 else
  869.                 {
  870.                   av=avma;p1=polx[vx];
  871.                   z=gsubst(x[l-1],v,y);
  872.                   for(i=l-1;i>3;i--)
  873.                     z=gadd(gmul(z,p1),gsubst(x[i-1],v,y));
  874.                   z=gmul(z,p1);p2=gsubst(x[2],v,y);
  875.                   tetpil=avma;
  876.                   z=gerepile(av,tetpil,gadd(p2,z));
  877.                 }
  878.               }
  879.             }
  880.           }
  881.         }
  882.       }
  883.       break;
  884.       
  885.     case 11:
  886.       ex=valp(x);vx=varn(x);
  887.       if(vx>v)
  888.       {
  889.         if(ty==19) z=gscalmat(x,ly-1);
  890.         else z=gcopy(x);
  891.         return z;
  892.       }
  893.       if(!signe(x))
  894.       {
  895.         if(vx<v) z=gcopy(x);
  896.         else
  897.         {
  898.           z=cgetg(3,11);
  899.           z[1]=0x8000+ex*gval(y,v);z[2]=zero;
  900.           setvarn(z,vx);
  901.         }
  902.         return z;
  903.       }
  904.       if(vx<v)
  905.       {
  906.         /* a ameliorer */
  907.         av1=avma;setvalp(x,0);p1=gconvsp(x);setvalp(x,ex);
  908.         p2=gsubst(p1,v,y);av2=avma;
  909.         z=tayl(p2,vx,lx-2);
  910.         if(ex)
  911.         {
  912.           p1=gpuigs(polx[vx],ex);
  913.           av2=avma;z=gerepile(av1,av2,gmul(z,p1));
  914.         }
  915.         else z=gerepile(av1,av2,z);
  916.         return z;
  917.       }
  918.       switch(ty)
  919.       {
  920.         case 11:
  921.     ey=valp(y);vy=varn(y);
  922.     if (ey<1)
  923.       {
  924.         z=cgetg(3,11);z[2]=zero;
  925.         z[1]=0x8000+ey*(ex+lx-2);
  926.         setvarn(z,vy);
  927.       }
  928.     else
  929.       {
  930.         l=(lx-2)*ey+2;
  931.         if (ex)
  932.           {if (l>ly) l=ly;}
  933.         else 
  934.           {
  935.         if (gcmp0(y)) l=ey+2;
  936.         else
  937.           {if (l>ly) l=ly+ey;}
  938.           }
  939.         if(vy!=vx)
  940.           {
  941.         av=avma;z=cgetg(2,11);z[1]=0x8000;setvarn(z,vy);
  942.         for(i=lx-1;i>=2;i--)
  943.           {p1=gmul(y,z);av2=avma;z=gadd(x[i],p1);}
  944.         if (ex)
  945.           {
  946.             p1=gpuigs(y,ex);av2=avma;
  947.             z=gerepile(av,av2,gmul(z,p1));
  948.           }
  949.         else z=gerepile(av,av2,z);
  950.           }
  951.         else
  952.           {
  953.         av1=avma;t=cgetg(ly,11);
  954.         av2=avma;z=cgetg(l,11);
  955.         z[2] =lcopy(x[2]);
  956.         for (i=3;i<l;i++) z[i]=zero;
  957.         for (i=2;i<ly;i++) t[i]=y[i];
  958.         
  959.         for (i=3,jb=ey;jb<=l-2;i++,jb+=ey)
  960.           {
  961.             for (j=jb+2;j<l;j++)
  962.               {
  963.             av4=avma;
  964.             p1=gmul(x[i],t[j-jb]);
  965.             av5=avma;
  966.             z[j]=lpile(av4,av5,ladd(z[j],p1));
  967.             if (j==jb+ey+1) av=avma;
  968.               }
  969.             for (j=l-1-jb-ey;j>1;j--)
  970.               {
  971.             av4=avma;p1=gzero;
  972.             for (k=2;k<j;k++)
  973.               p1=gadd(p1,gmul(t[j-k+2],y[k]));
  974.             p2=gmul(t[2],y[j]);
  975.             av5=avma;
  976.             t[j]=lpile(av4,av5,gadd(p1,p2));
  977.               }
  978.             if (i>3)
  979.               {
  980.             av7=avma;decal=lpile(av3,av6,0);
  981.             for (k=jb+2;k<l;k++)
  982.               if (adecaler(z[k],av6,av7)) z[k]+=decal;
  983.             for (k=2;k<l-jb-ey+2;k++)
  984.               if (adecaler(t[k],av6,av7)) t[k]+=decal;
  985.             av3=av+decal;
  986.               }
  987.             else av3=av;
  988.             av6=avma;
  989.           }
  990.         z[1]=0x1008000;setvarn(z,varn(y));
  991.         if (ex)
  992.           {
  993.             if (l<ly) setlg(y,l);
  994.             p1=gpuigs(y,ex);
  995.             av2=avma;
  996.             z=gerepile(av1,av2,gmul(z,p1));
  997.             if (l<ly) setlg(y,ly);
  998.           }
  999.         else z=gerepile(av1,av2,z);
  1000.           }
  1001.       }
  1002.           break;
  1003.           
  1004.         case 10:
  1005.         case 13:
  1006.         case 14:
  1007.           if(gcmp0(y))
  1008.           {
  1009.             z=cgetg(lx,tx);z[1]=0x1008000;
  1010.             z[2]=lcopy(x[2]);setvarn(z,v);
  1011.             for(i=3;i<lx;i++) z[i]=zero;
  1012.           }
  1013.           else
  1014.           {
  1015.             e=gval(y,vy=gvar(y));if(e<=0) err(gsubser5);
  1016.         av=avma;p1=gconvsp(x);p2=gsubst(p1,v,y);tetpil=avma;
  1017.         z=gerepile(av,tetpil,tayl(p2,vy,e*(lx-2+ex)));
  1018.           }
  1019.           break;
  1020.         default: err(gsubser4);
  1021.         }
  1022.       break;
  1023.     case 9:
  1024.       z=cgetg(lx,tx);z[1]=lsubst(x[1],v,y);av=avma;
  1025.       p1=gsubst(x[2],v,y);tetpil=avma;
  1026.       z[2]=lpile(av,tetpil,gmod(p1,z[1]));
  1027.       break;
  1028. /*  case 9: err(gsubser1); */
  1029.     case 13:
  1030.     case 14:
  1031.       av=avma;p1=gsubst(x[1],v,y);p2=gsubst(x[2],v,y);
  1032.       tetpil=avma;z=gerepile(av,tetpil,gdiv(p1,p2));
  1033.       break;
  1034.     case 17:
  1035.     case 18:
  1036.     case 19:
  1037.       z=cgetg(lx,tx);
  1038.       for(i=1;i<lx;i++)
  1039.       z[i]=lsubst(x[i],v,y);
  1040.   }
  1041.   return z;
  1042. }
  1043.  
  1044. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1045. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1046. /*                                                                 */
  1047. /*                SERIE RECIPROQUE D'UNE SERIE                     */
  1048. /*                                                                 */
  1049. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1050. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1051.  
  1052. GEN     recip(x)
  1053.      
  1054.      GEN     x;
  1055.      
  1056. {
  1057.   
  1058.   long  lx,av1,av2,av3,av4,av5,av6,av7,av8,i,j,k,v,decal;
  1059.   GEN   a,y,u,p1,p2;
  1060.   
  1061.   if ((typ(x)!=11) || (valp(x)!=1)) err(reciper);
  1062.   /* attention au coeff directeur */
  1063.   a=(GEN)x[2];v=varn(x);
  1064.   if (gcmp1(a))
  1065.   {
  1066.     lx=lg(x);av1=avma;u=cgetg(lx,11);
  1067.     av2=avma;y=cgetg(lx,11);
  1068.     u[2]=un;y[2]=un;
  1069.     av5=avma;u[3]=lmulsg(-2,x[3]);
  1070.     av6=avma;y[3]=lneg(x[3]);
  1071.     av7=avma;i=2;
  1072.     while (++i<lx-1)
  1073.     {
  1074.       for (j=3;j<i+1;j++)
  1075.       {
  1076.         av3=avma;p1=(GEN)u[j];
  1077.         for (k=j-1;k>2;k--)
  1078.           p1=gsub(p1,gmul(u[k],x[j-k+2]));
  1079.         av4=avma;
  1080.         u[j]=lpile(av3,av4,lsub(p1,x[j]));
  1081.       }
  1082.       av3=avma;p1=gmulsg(i,x[i+1]);
  1083.       for (k=2;k<i;k++)
  1084.       {
  1085.         p2=gmul(x[k+1],u[i-k+2]);
  1086.         p1=gadd(p1,gmulsg(k,p2) );
  1087.       }
  1088.       av4=avma;u[i+1]=lpile(av3,av4,lneg(p1));
  1089.       av8=avma;decal=lpile(av5,av6,0);
  1090.       for (k=3;k<i+2;k++)
  1091.         if (adecaler(u[k],av6,av8)) u[k]+=decal;
  1092.       if(adecaler(y[i],av6,av8)) y[i]+=decal;
  1093.       av6=avma;av5=av7+decal;
  1094.       y[i+1]=ldivgs(u[i+1],i);av7=avma;
  1095.     }
  1096.     y[lx-1]=lpile(av5,av6,y[lx-1]);
  1097.     y[1]=0x1008001;y=gerepile(av1,av2,y);
  1098.     setvarn(y,v);
  1099.   }
  1100.   else
  1101.   {
  1102.     p1=(GEN)x[2];av1=avma;y=gdiv(x,p1);
  1103.     y=recip(y);p2=gdiv(polx[v],p1);av2=avma;
  1104.     y=gerepile(av1,av2,gsubst(y,v,p2));
  1105.   }
  1106.   return y;
  1107. }
  1108.  
  1109.  
  1110. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1111. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1112. /*                                                                 */
  1113. /*                    DERIVATION ET INTEGRATION                    */
  1114. /*                                                                 */
  1115. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1116. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1117.  
  1118. GEN     deriv(x,v)
  1119.      
  1120.      GEN     x;
  1121.      long    v;
  1122.      
  1123. {
  1124.   long  lx,ly,vx,tx,e,i,j,tetpil,l,l1,f;
  1125.   GEN   y,p1,p2;
  1126.   
  1127.   tx=typ(x);if(tx<9) return gzero;
  1128.   switch(tx)
  1129.     {
  1130.     case 9: if(v<=varn(x[1])) return gzero;
  1131.       y=cgetg(3,9);y[1]=copyifstack(x[1]);y[2]=lderiv(x[2],v);
  1132.       return y;break;
  1133.     case 10: vx=varn(x);
  1134.       if(vx>v) return gzero;
  1135.       lx=lgef(x)-1;
  1136.       if(vx<v)
  1137.       {
  1138.         l=avma;
  1139.         for(i=lx;(i>=2)&&(isexactzero(p1=deriv(x[i],v)));avma=l,i--);
  1140.         y=cgetg(i+1,10);
  1141.         if(i==1) y[1]=2;
  1142.         else
  1143.         {
  1144.           y[1]=0x1000001+i;y[i]=(long)p1;f=gcmp0(p1);
  1145.           for(j=2;j<i;j++) {y[j]=lderiv(x[j],v);if(f) f=gcmp0(y[j]);}
  1146.       if(f) setsigne(y,0);
  1147.         }
  1148.         setvarn(y,vx);
  1149.         return y;
  1150.       }
  1151.       if(lx<3) y=gzero;
  1152.       else
  1153.       {
  1154.         l=avma;
  1155.         for(i=lx-1;(i>=2)&&isexactzero(p1=gmulsg(i-1,x[i+1]));avma=l,i--);
  1156.         y=cgetg(i+1,tx);
  1157.         if(i==1) y[1]=2;
  1158.         else
  1159.         {
  1160.           y[1]=0x1000001+i;y[i]=(long)p1;f=gcmp0(p1);
  1161.           for(i--;i>=2;i--) {y[i]=lmulsg(i-1,x[i+1]);if(f) f=gcmp0(y[i]);}
  1162.       if(f) setsigne(y,0);
  1163.         }
  1164.         setvarn(y,v);
  1165.       }
  1166.       break;
  1167.     case 11: vx=varn(x);
  1168.       if(vx>v) return gzero;
  1169.       lx=lg(x);e=valp(x);
  1170.       if(vx<v)
  1171.       {
  1172.         if(!signe(x)) y=gcopy(x);
  1173.         else
  1174.         {
  1175.           l=avma;
  1176.           for(i=2;(i<lx)&&(gcmp0(p1=deriv(x[i],v)));avma=l,i++);
  1177.           if(i==lx) y=ggrando(polx[vx],e+lx-2);
  1178.           else
  1179.           {
  1180.             y=cgetg(lx-i+2,11);y[1]=0x7ffe +e+i;
  1181.             setvarn(y,vx);y[2]=(long)p1;
  1182.             for(j=3;j<=lx-i+1;j++)
  1183.               y[j]=lderiv(x[i+j-2],v);
  1184.           }
  1185.         }
  1186.         return y;
  1187.       }
  1188.       ly=lx-1;if(ly<3) ly=3;
  1189.       if(gcmp0(x))
  1190.         {y=cgetg(3,tx);y[1]=0x7fff+e;}
  1191.       else
  1192.       {
  1193.         if(e)
  1194.         {
  1195.           y=cgetg(lx,tx);
  1196.           for(i=2;i<lx;i++)
  1197.             y[i]=lmulsg(i+e-2,x[i]);
  1198.           y[1]=0x1008000+e-1;
  1199.         }
  1200.         else
  1201.         {
  1202.           i=3;
  1203.           while((i<lx)&&gcmp0(x[i])) i++;
  1204.           if(i==lx)
  1205.             {y=cgetg(ly,tx);y[1]=0x7ffd+lx;}
  1206.           else
  1207.           {
  1208.             ly=ly-i+3;y=cgetg(ly,tx);
  1209.             y[1]=0x1007ffd+i;
  1210.             for(j=2;j<ly;j++)
  1211.               y[j]=lmulsg(j+i-4,x[i+j-2]);
  1212.           }
  1213.         }
  1214.       }
  1215.       setvarn(y,vx);
  1216.       break;
  1217.     case 13:
  1218.     case 14: l=avma;y=cgetg(3,tx);y[2]=lmul(x[2],x[2]);
  1219.       l1=avma;p1=gmul(x[2],deriv(x[1],v));p2=gmul(x[1],deriv(x[2],v));
  1220.       tetpil=avma;p1=gsub(p1,p2);y[1]=lpile(l1,tetpil,p1);
  1221.       if(tx==13) {tetpil=avma;y=gerepile(l,tetpil,gred(y));}
  1222.       break;
  1223.     case 17:
  1224.     case 18:
  1225.     case 19: lx=lg(x);y=cgetg(lx,tx);
  1226.       for(i=1;i<lx;i++)
  1227.         y[i]=lderiv(x[i],v);
  1228.       break;
  1229.     default: break;
  1230.     }
  1231.   return y;
  1232. }
  1233.  
  1234.  
  1235. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1236. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1237. /*                                                                 */
  1238. /*                    INTEGRATION FORMELLE                         */
  1239. /*                                                                 */
  1240. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1241. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1242.  
  1243. GEN     integ(x,v)
  1244.      
  1245.      GEN     x;
  1246.      long    v;
  1247.      
  1248. {
  1249.   long  lx,tx,e,i,j,vx;
  1250.   GEN   y;
  1251.   
  1252.   tx=typ(x);
  1253.   if((tx<9)||((tx==9)&&(v<=varn(x[1]))))
  1254.   {
  1255.     if(gcmp0(x)) y=gzero;
  1256.     else
  1257.     {
  1258.       y=cgetg(4,10);y[1]=0x1000004;setvarn(y,v);
  1259.       y[2]=zero;y[3]=lcopy(x);
  1260.     }
  1261.     return y;
  1262.   }
  1263.   switch(tx)
  1264.     {
  1265.     case 9: y=cgetg(3,9);y[1]=copyifstack(x[1]);y[2]=linteg(x[2],v);
  1266.       return y;break;
  1267.     case 10: vx=varn(x);lx=lgef(x);
  1268.       if(lx==2)
  1269.       {
  1270.         y=cgetg(3,10);y[1]=2;
  1271.         if(vx>v) setvarn(y,v);else setvarn(y,vx);
  1272.         return y;
  1273.       }
  1274.       if(vx>v)
  1275.       {
  1276.         y=cgetg(4,10);y[1]=(signe(x))?0x1000004:4;setvarn(y,v);
  1277.         y[2]=zero;y[3]=lcopy(x);
  1278.         return y;
  1279.       }
  1280.       if(vx<v)
  1281.       {
  1282.         y=cgetg(lx,10);y[1]=x[1];
  1283.         for(i=2;i<lx;i++)
  1284.           y[i]=linteg(x[i],v);
  1285.         return y;
  1286.       }
  1287.       y=cgetg(lx+1,tx);y[2]=zero;
  1288.       for(i=3;i<=lx;i++)
  1289.         y[i]=ldivgs(x[i-1],i-2);
  1290.       y[1]=signe(x)?0x1000001+lx:1+lx;setvarn(y,v);
  1291.       break;
  1292.     case 11: lx=lg(x);e=valp(x);vx=varn(x);
  1293.       if(!signe(x))
  1294.       {
  1295.         y=cgetg(3,tx);
  1296.         if(vx==v) y[1]=0x8001+e;
  1297.         else y[1]=x[1];
  1298.         if(vx>v) setvarn(y,v);
  1299.         return y;
  1300.       }
  1301.       if(vx>v)
  1302.       {
  1303.         y=cgetg(4,10);y[1]=0x1000004;setvarn(y,v);
  1304.         y[2]=zero;y[3]=lcopy(x);
  1305.         return y;
  1306.       }
  1307.       if(vx<v)
  1308.       {
  1309.         y=cgetg(lx,tx);y[1]=x[1];
  1310.         for(i=2;i<lx;i++)
  1311.           y[i]=linteg(x[i],v);
  1312.         return y;
  1313.       }
  1314.       y=cgetg(lx,tx);
  1315.       for(i=2;i<lx;i++)
  1316.       {
  1317.         if(!(j=i+e-1)) err(inter2);
  1318.         y[i]=ldivgs(x[i],j);
  1319.       }
  1320.       y[1]=(x[1])+1;
  1321.       break;
  1322.     case 13:
  1323.     case 14: err(impl,"integration of a rational function");
  1324.     case 17:
  1325.     case 18:
  1326.     case 19: lx=lg(x);y=cgetg(lx,tx);
  1327.       for(i=1;i<lx;i++)
  1328.         y[i]=linteg(x[i],v);
  1329.       break;
  1330.     default: err(inter1);
  1331.   }
  1332.   return y;
  1333. }
  1334.  
  1335. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1336. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1337. /*                                                                 */
  1338. /*                    PARTIES ENTIERES                             */
  1339. /*                                                                 */
  1340. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1341. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1342.  
  1343. GEN     gfloor(x)
  1344.      
  1345.      GEN     x;
  1346.      
  1347. {
  1348.   GEN     y,p1;
  1349.   long  i,lx,tx,av,tetpil;
  1350.   
  1351.   switch(tx=typ(x))
  1352.   {
  1353.     case 1 :
  1354.     case 10: y=gcopy(x);break;
  1355.     case 2 : y=mpent(x);break;
  1356.     case 4 :
  1357.     case 5 : av=avma;y=dvmdii(x[1],x[2],&p1);
  1358.       i=!gcmp0(p1);cgiv(p1);
  1359.       if(i&&(gsigne(x)<0))
  1360.       {
  1361.         tetpil=avma;y=gerepile(av,tetpil,gsubgs(y,1));
  1362.       }
  1363.       break;
  1364.     case 13:
  1365.     case 14: y=gdeuc(x[1],x[2]);
  1366.       break;
  1367.     case 17:
  1368.     case 18:
  1369.     case 19: lx=lg(x);y=cgetg(lx,tx);
  1370.       for(i=1;i<lx;i++)
  1371.         y[i]=lfloor(x[i]);
  1372.       break;
  1373.     default: err(flooer);
  1374.   }
  1375.   return y;
  1376. }
  1377.  
  1378. GEN     gfrac(x)
  1379.      
  1380.      GEN        x;
  1381.      
  1382. {
  1383.   long  av,tetpil;
  1384.   GEN   p1;
  1385.   
  1386.   av=avma;p1=gfloor(x);tetpil=avma;
  1387.   return gerepile(av,tetpil,gsub(x,p1));
  1388. }
  1389.  
  1390. GEN     gceil(x)
  1391.      
  1392.      GEN     x;
  1393.      
  1394. {
  1395.   GEN     y,p1;
  1396.   long  i,lx,tx=typ(x),av,tetpil;
  1397.   
  1398.   switch(tx)
  1399.   {
  1400.     case 1 :
  1401.     case 10: y=gcopy(x);break;
  1402.     case 2 : av=avma;y=mpent(x);
  1403.       if(!gegal(x,y))
  1404.       {
  1405.         tetpil=avma;y=gerepile(av,tetpil,gaddsg(1,y));
  1406.       }
  1407.       break;
  1408.     case 4 :
  1409.     case 5 : av=avma;y=dvmdii(x[1],x[2],&p1);
  1410.       i=!gcmp0(p1);cgiv(p1);
  1411.       if(i&&(gsigne(x)>0))
  1412.       {
  1413.         tetpil=avma;y=gerepile(av,tetpil,gaddsg(1,y));
  1414.       }
  1415.       break;
  1416.     case 13:
  1417.     case 14: y=gdeuc(x[1],x[2]);
  1418.       break;
  1419.     case 17:
  1420.     case 18:
  1421.     case 19: lx=lg(x);y=cgetg(lx,tx);
  1422.       for(i=1;i<lx;i++)
  1423.         y[i]=lceil(x[i]);
  1424.       break;
  1425.     default: err(ceiler);
  1426.   }
  1427.   return y;
  1428. }
  1429.  
  1430. GEN     ground(x)
  1431.      
  1432.      GEN     x;
  1433.      
  1434. {
  1435.   GEN     y,p1;
  1436.   long  i,lx=lg(x),tx=typ(x),av,tetpil;
  1437.   
  1438.   switch(tx)
  1439.   {
  1440.     case 1 :
  1441.     case 3 :
  1442.     case 8 : y=gcopy(x);break;
  1443.     case 2 : 
  1444.     case 4 :
  1445.     case 5 : av=avma;p1=gadd(ghalf,x);tetpil=avma;
  1446.       y=gerepile(av,tetpil,gfloor(p1));break;
  1447.     case 9 : y=cgetg(3,9);y[1]=copyifstack(x[1]);y[2]=lround(x[2]);
  1448.       break;
  1449.     case 10: lx=lgef(x);
  1450.     case 6 :
  1451.     case 11:
  1452.     case 13:
  1453.     case 14: 
  1454.     case 17:
  1455.     case 18:
  1456.     case 19: lx=lg(x);y=cgetg(lx,tx);
  1457.       for(i=1;i<lontyp[tx];i++) y[i]=x[i];
  1458.       for(i=lontyp[tx];i<lx;i++) y[i]=lround(x[i]);
  1459.       break;
  1460.     default: err(rounder);
  1461.   }
  1462.   return y;
  1463. }
  1464.  
  1465. GEN     grndtoi(x,e)
  1466.      
  1467.      GEN     x;
  1468.      long   *e;
  1469.      
  1470. {
  1471.   GEN     y,p1;
  1472.   long  i,lx=lg(x),tx=typ(x),av,tetpil,ex,e1;
  1473.   
  1474.   *e= -1048576;
  1475.   switch(tx)
  1476.   {
  1477.     case 1 :
  1478.     case 3 :
  1479.     case 8 : y=gcopy(x);break;
  1480.     case 2 : av=avma;p1=gadd(ghalf,x);
  1481.       ex=expo(p1);*e=ex-(lx-2)*32;
  1482.       if(ex<0)
  1483.       {
  1484.     if(signe(p1)>=0) {avma=av;y=gzero;}
  1485.     else {tetpil=avma;y=gerepile(av,tetpil,gneg(un));}
  1486.       }
  1487.       else
  1488.       {
  1489.         settyp(p1,1);setlgef(p1,lx);
  1490.         if(signe(x)>=0)
  1491.           {tetpil=avma;y=gerepile(av,tetpil,shifti(p1,*e+1));}
  1492.         else
  1493.         {
  1494.           y=shifti(p1,*e+1);tetpil=avma;
  1495.           y=gerepile(av,tetpil,gaddgs(y,-1));
  1496.         }
  1497.       }
  1498.       break;
  1499.     case 4 :
  1500.     case 5 : av=avma;p1=gadd(ghalf,x);tetpil=avma;
  1501.       y=gerepile(av,tetpil,gfloor(p1));break;
  1502.     case 9 : y=cgetg(3,9);y[1]=copyifstack(x[1]);y[2]=lrndtoi(x[2],e);
  1503.       break;
  1504.     case 10: lx=lgef(x);
  1505.     case 6 :
  1506.     case 11:
  1507.     case 13:
  1508.     case 14: 
  1509.     case 17:
  1510.     case 18:
  1511.     case 19: lx=lg(x);y=cgetg(lx,tx);
  1512.       for(i=1;i<lontyp[tx];i++)
  1513.         y[i]=x[i];
  1514.       for(i=lontyp[tx];i<lx;i++)
  1515.       {
  1516.         y[i]=lrndtoi(x[i],&e1);
  1517.         if(e1>*e) *e=e1;
  1518.       }
  1519.       break;
  1520.     default: err(rndtoier);
  1521.   }
  1522.   return y;
  1523. }
  1524.  
  1525. long rounderror(x)
  1526.      GEN x;
  1527.      
  1528. {
  1529.   long e,av=avma;
  1530.  
  1531.   grndtoi(x,&e);avma=av;
  1532.   e*=L2SL10;return e;
  1533. }
  1534.  
  1535. GEN     gcvtoi(x,e)
  1536.      
  1537.      GEN     x;
  1538.      long    *e;
  1539.      
  1540.      /* la variable formelle e,representant le nombre de bits d'erreur */
  1541.      /* sur la partie entiere,n'est utilisee que dans le cas 2 (reel) */
  1542.      
  1543. {
  1544.   GEN     y,p1;
  1545.   long  tx=typ(x),lx=lg(x),i,ex,av,tetpil,e1;
  1546.   
  1547.   *e= -1048576;
  1548.   switch(tx)
  1549.   {
  1550.     case 1 :
  1551.     case 10: y=gcopy(x);break;
  1552.     case 2 : ex=expo(x);*e=ex-(lx-2)*32;
  1553.       if(ex<0) y=gzero;
  1554.       else
  1555.       {
  1556.         av=avma;p1=gcopy(x);settyp(p1,1);setlgef(p1,lx);
  1557.         tetpil=avma;y=gerepile(av,tetpil,shifti(p1,*e+1));
  1558.       }
  1559.       break;
  1560.     case 4 :
  1561.     case 5 : y=divii(x[1],x[2]);  break;
  1562.     case 7 : y=gconvpe(x);break;
  1563.     case 13:
  1564.     case 14: y=gdeuc(x[1],x[2]);  break;
  1565.     case 11: y=gconvsp(x);break;
  1566.     case 17:
  1567.     case 18:
  1568.     case 19: y=cgetg(lx,tx);
  1569.       for(i=1;i<lx;i++)
  1570.       {
  1571.         y[i]=lcvtoi(x[i],&e1);
  1572.         if(e1>*e) *e=e1;
  1573.       }
  1574.       break;
  1575.     default: err(cvtoier);
  1576.   }
  1577.   return y;
  1578. }
  1579.  
  1580. GEN     gtrunc(x)
  1581.      
  1582.      GEN     x;
  1583.      
  1584. {
  1585.   GEN     y;
  1586.   long  tx=typ(x),lx=lg(x),i;
  1587.   switch(tx)
  1588.   {
  1589.     case 1 :
  1590.     case 10: y=gcopy(x);break;
  1591.     case 2 : y=mptrunc(x);break;
  1592.     case 4 :
  1593.     case 5 : y=divii(x[1],x[2]); break;
  1594.     case 7 : y=gconvpe(x);break;
  1595.     case 13:
  1596.     case 14: y=gdeuc(x[1],x[2]); break;
  1597.     case 11: y=gconvsp(x);break;
  1598.     case 17:
  1599.     case 18:
  1600.     case 19: y=cgetg(lx,tx);
  1601.       for(i=1;i<lx;i++)
  1602.         y[i]=ltrunc(x[i]);
  1603.       break;
  1604.     default: err(truncer);
  1605.   }
  1606.   return y;
  1607. }
  1608.  
  1609. GEN gtopoly(x,v)
  1610.      GEN x;
  1611.      long v;
  1612. {
  1613.   long tx=typ(x),lx,i,j;
  1614.   GEN y;
  1615.   
  1616.   if(gcmp0(x)) {y=cgetg(2,10);y[1]=2;setvarn(y,v);return y;}
  1617.   if(tx<10)
  1618.     {
  1619.       y=cgetg(3,10);y[1]=0x1000003;y[2]=lcopy(x);
  1620.       setvarn(y,v);return y;
  1621.     }
  1622.   switch(tx)
  1623.     {
  1624.     case 10: y=gcopy(x);break;
  1625.     case 11: y=gconvsp(x);break;
  1626.     case 13:
  1627.     case 14: y=gdeuc(x[1],x[2]);break;
  1628.     case 15:
  1629.     case 16:
  1630.     case 17:
  1631.     case 18:
  1632.     case 19: lx=lg(x);i=1;while((i<lx)&&gcmp0(x[i])) i++;
  1633.       y=cgetg(lx-i+2,10);y[1]=0x1000002+lx-i;
  1634.       for(j=2;j<=lx-i+1;j++) y[j]=lcopy(x[lx+1-j]);
  1635.     }
  1636.   setvarn(y,v);return y;
  1637. }
  1638.  
  1639. GEN gtopolyrev(x,v)
  1640.      GEN x;
  1641.      long v;
  1642. {
  1643.   long tx=typ(x),lx,i,j;
  1644.   GEN y;
  1645.   
  1646.   if(gcmp0(x)) {y=cgetg(2,10);y[1]=2;setvarn(y,v);return y;}
  1647.   if(tx<10)
  1648.     {
  1649.       y=cgetg(3,10);y[1]=0x1000003;y[2]=lcopy(x);
  1650.       setvarn(y,v);return y;
  1651.     }
  1652.   switch(tx)
  1653.     {
  1654.     case 10: y=gcopy(x);break;
  1655.     case 11: y=gconvsp(x);break;
  1656.     case 13:
  1657.     case 14: y=gdeuc(x[1],x[2]);break;
  1658.     case 15:
  1659.     case 16:
  1660.     case 17:
  1661.     case 18:
  1662.     case 19: lx=lg(x);i=1;while((i<lx)&&gcmp0(x[lx-i])) i++;
  1663.       y=cgetg(lx-i+2,10);y[1]=0x1000002+lx-i;
  1664.       for(j=2;j<=lx-i+1;j++) y[j]=lcopy(x[j-1]);
  1665.     }
  1666.   setvarn(y,v);return y;
  1667. }
  1668.  
  1669. GEN gtoser(x,v)
  1670.      GEN x;
  1671.      long v;
  1672. {
  1673.   long tx=typ(x),lx,i,j,l,av,tetpil;
  1674.   GEN y,p1,p2;
  1675.   
  1676.   if(tx==11) {y=gcopy(x);setvarn(y,v);return y;}
  1677.   if(gcmp0(x)) {y=cgetg(2,11);y[1]=0x8000+precdl;setvarn(y,v);return y;}
  1678.   if(tx<10)
  1679.     {
  1680.       y=cgetg(precdl+2,11);y[1]=0x1008000;y[2]=lcopy(x);
  1681.       for(i=3;i<=precdl+1;i++) y[i]=zero;
  1682.       setvarn(y,v);return y;
  1683.     }
  1684.   switch(tx)
  1685.     {
  1686.     case 10: lx=lgef(x);i=2;while((i<lx)&&gcmp0(x[i])) i++;
  1687.       l=lx-i;if(precdl>l) l=precdl;
  1688.       y=cgetg(l+2,11);y[1]=0x1008000+i-2;
  1689.       for(j=2;j<=lx-i+1;j++) y[j]=lcopy(x[j+i-2]);
  1690.       for(j=lx-i+2;j<=l+1;j++) y[j]=zero;
  1691.       break;
  1692.     case 13:
  1693.     case 14: av=avma;p1=gtoser(x[1],v);p2=gtoser(x[2],v);
  1694.       tetpil=avma;y=gerepile(av,tetpil,gdiv(p1,p2));break;
  1695.     case 15:
  1696.     case 16:
  1697.     case 17:
  1698.     case 18:
  1699.     case 19: lx=lg(x);i=1;while((i<lx)&&gcmp0(x[i])) i++;
  1700.       y=cgetg(lx-i+2,11);y[1]=0x1008000+i-1;
  1701.       for(j=2;j<=lx-i+1;j++) y[j]=lcopy(x[j+i-2]);
  1702.     }
  1703.   setvarn(y,v);return y;
  1704. }
  1705.  
  1706. GEN gtovec(x)
  1707.      GEN x;
  1708. {
  1709.   long tx=typ(x),lx,i;
  1710.   GEN y;
  1711.   
  1712.   if((tx<10)||(tx==13)||(tx==14))
  1713.     {y=cgetg(2,17);y[1]=lcopy(x);return y;}
  1714.   if(tx>=15)
  1715.     {
  1716.       lx=lg(x);y=cgetg(lx,17);
  1717.       for(i=1;i<lx;i++) y[i]=lcopy(x[i]);
  1718.       return y;
  1719.     }
  1720.   if(tx==10)
  1721.     {
  1722.       lx=lgef(x);y=cgetg(lx-1,17);
  1723.       for(i=1;i<=lx-2;i++) y[i]=lcopy(x[lx-i]);
  1724.       return y;
  1725.     }
  1726.   if(!signe(x)) {y=cgetg(2,17);y[1]=zero;return y;}
  1727.   lx=lg(x);y=cgetg(lx-1,17);
  1728.   for(i=1;i<=lx-2;i++) y[i]=lcopy(x[i+1]);
  1729.   return y;
  1730. }
  1731.       
  1732. GEN     compo(x,n)
  1733.      
  1734.      GEN     x;
  1735.      long    n;
  1736.      
  1737. {
  1738.   long l,lx=lg(x),tx=typ(x);
  1739.  
  1740.   if((tx==10)&&((n+1)>=lgef(x))) return gzero;
  1741.   if((tx==11)&&(!signe(x))) return gzero;
  1742.   l=lontyp[tx]+n-1;
  1743.   if((l>=lx) || (n<1)) err(compoer1);
  1744.   return gcopy(x[l]);
  1745. }
  1746.  
  1747. GEN truecoeff(x,n)
  1748.      GEN x;
  1749.      long n;
  1750. {
  1751.   long tx=typ(x),lx,ex;
  1752.  
  1753.   if(tx<10)
  1754.     {if(n) return gzero;else return gcopy(x);}
  1755.   switch(tx)
  1756.     {
  1757.     case 15: case 16: case 17: case 18: case 19:
  1758.       if((n<=0)||(n>=lg(x))) err(compoer1);
  1759.       return gcopy(x[n]);break;
  1760.     case 10:
  1761.       if((n<0)||(n>=lgef(x)-2)) return gzero;
  1762.       return gcopy(x[n+2]);break;
  1763.     case 11:
  1764.       if(!signe(x)) return gzero;
  1765.       lx=lg(x);ex=valp(x);
  1766.       if(n<ex) return gzero;
  1767.       if(n>=ex+lx-2) err(compoer1);
  1768.       return gcopy(x[n-ex+2]);break;
  1769.     default: err(compoer1);
  1770.     }
  1771. }
  1772.  
  1773. GEN  denom(x)
  1774.      GEN x;
  1775. {
  1776.   long i,av,tetpil;
  1777.   GEN s,t;
  1778.  
  1779.   switch(typ(x))
  1780.     {
  1781.     case 1: case 2: case 3: case 7: case 11: return gun;
  1782.     case 4: case 5: return absi(x[2]);
  1783.     case 6: av=avma;t=denom(x[1]);s=denom(x[2]);tetpil=avma;
  1784.       return gerepile(av,tetpil,glcm(s,t));
  1785.     case 8: av=avma;t=denom(x[2]);s=denom(x[3]);tetpil=avma;
  1786.       return gerepile(av,tetpil,glcm(s,t));
  1787.     case 9: return denom(x[2]);
  1788.     case 13: case 14: return gcopy(x[2]);
  1789.     case 10: return polun[varn(x)];
  1790.     case 17: case 18: case 19: av=tetpil=avma;
  1791.       s=denom(x[1]);
  1792.       for(i=2;(i<lg(x));i++) {t=denom(x[i]);tetpil=avma;s=glcm(s,t);}
  1793.       return gerepile(av,tetpil,s);
  1794.     default: err(denomer1);
  1795.     }
  1796. }
  1797.  
  1798. GEN  numer(x)
  1799.      GEN x;
  1800. {
  1801.   long av,tetpil;
  1802.   GEN s;
  1803.  
  1804.   switch(typ(x))
  1805.     {
  1806.     case 1: case 10: case 2: case 3: case 7: case 11: return gcopy(x);
  1807.     case 4: case 5: return (signe(x[2])>0) ? gcopy(x[1]) : gneg(x[1]);
  1808.     case 9: av=avma;s=numer(x[2]);tetpil=avma;return gerepile(av,tetpil,gmodulcp(s,x[1]));
  1809.     case 13: case 14: return gcopy(x[1]);
  1810.     case 6: case 8: case 17: case 18: case 19: av=avma;s=denom(x);tetpil=avma;
  1811.       return gerepile(av,tetpil,gmul(s,x));
  1812.     default: err(numerer1);
  1813.     }
  1814. }
  1815.  
  1816. GEN lift(x)
  1817.      GEN x;
  1818. {
  1819.   long lx,tx=typ(x),i;
  1820.   GEN y;
  1821.  
  1822.   switch(tx)
  1823.     {
  1824.     case 1: return gcopy(x);
  1825.     case 3:
  1826.     case 9: return gcopy(x[2]);
  1827.     case 11: if(!signe(x)) return gcopy(x);
  1828.       y=cgetg(lx=lg(x),tx);y[1]=x[1];
  1829.       for(i=2;i<lx;i++) y[i]=(long)lift(x[i]);break;
  1830.     case 4:
  1831.     case 5:
  1832.     case 6:
  1833.     case 13:
  1834.     case 14:
  1835.     case 17:
  1836.     case 18:
  1837.     case 19: y=cgetg(lx=lg(x),tx);
  1838.       for(i=1;i<lx;i++) y[i]=(long)lift(x[i]);
  1839.       break;
  1840.     case 10: y=cgetg(lx=lgef(x),tx);y[1]=x[1];
  1841.       for(i=2;i<lx;i++) y[i]=(long)lift(x[i]);break;
  1842.     case 8: y=cgetg(4,tx);y[1]=copyifstack(x[1]);
  1843.       for(i=2;i<4;i++) y[i]=(long)lift(x[i]);break;
  1844.     default: err(lifter1);
  1845.     }
  1846.   return y;
  1847. }
  1848.  
  1849. GEN centerlift(x)
  1850.      GEN x;
  1851. {
  1852.   long lx,tx=typ(x),i,av;
  1853.   GEN y;
  1854.  
  1855.   switch(tx)
  1856.     {
  1857.     case 1: return gcopy(x);
  1858.     case 3: av=avma;i=cmpii(shifti(x[2],1),x[1]);avma=av;
  1859.       y=(i>0)?subii(x[2],x[1]):gcopy(x[2]);break;
  1860.     case 9: y=gcopy(x[2]);break;
  1861.     case 11: if(!signe(x)) return gcopy(x);
  1862.       y=cgetg(lx=lg(x),tx);y[1]=x[1];
  1863.       for(i=2;i<lx;i++) y[i]=(long)centerlift(x[i]);break;
  1864.     case 4:
  1865.     case 5:
  1866.     case 6:
  1867.     case 13:
  1868.     case 14:
  1869.     case 17:
  1870.     case 18:
  1871.     case 19: y=cgetg(lx=lg(x),tx);
  1872.       for(i=1;i<lx;i++) y[i]=(long)centerlift(x[i]);
  1873.       break;
  1874.     case 10: y=cgetg(lx=lgef(x),tx);y[1]=x[1];
  1875.       for(i=2;i<lx;i++) y[i]=(long)centerlift(x[i]);break;
  1876.     case 8: y=cgetg(4,tx);y[1]=copyifstack(x[1]);
  1877.       for(i=2;i<4;i++) y[i]=(long)centerlift(x[i]);break;
  1878.     default: err(lifter1);
  1879.     }
  1880.   return y;
  1881. }
  1882.  
  1883. GEN glt(x,y) GEN x,y; { return gcmp(x,y)<0 ? gun : gzero;}
  1884. GEN gle(x,y) GEN x,y; { return gcmp(x,y)<=0 ? gun : gzero;}
  1885. GEN gge(x,y) GEN x,y; { return gcmp(x,y)>=0 ? gun : gzero;}
  1886. GEN ggt(x,y) GEN x,y; { return gcmp(x,y)>0 ? gun : gzero;}
  1887. GEN geq(x,y) GEN x,y; { return gegal(x,y) ? gun : gzero;}
  1888. GEN gne(x,y) GEN x,y; { return gegal(x,y) ? gzero : gun;}
  1889. GEN gand(x,y) GEN x,y; { return gcmp0(x)||gcmp0(y) ? gzero : gun;}
  1890. GEN gor(x,y) GEN x,y; { return gcmp0(x)&&gcmp0(y) ? gzero : gun;}
  1891.  
  1892. GEN glength(x)
  1893.   GEN x;
  1894. {
  1895.   switch(typ(x))
  1896.   {
  1897.     case 1:
  1898.     case 10: return stoi(lgef(x)-2);
  1899.     case 2: return stoi(lg(x)-2);
  1900.     case 11: return signe(x) ? stoi(lg(x)-2) : gzero;
  1901.     default: return stoi(lg(x)-lontyp[typ(x)]);
  1902.   }
  1903. }
  1904.  
  1905. GEN matsize(x)
  1906.      GEN x;
  1907. {
  1908.   GEN y;
  1909.  
  1910.   switch(typ(x))
  1911.     {
  1912.     case 17: y=cgetg(3,17);y[1]=un;y[2]=lstoi(lg(x)-1);break;
  1913.     case 18: y=cgetg(3,17);y[1]=lstoi(lg(x)-1);y[2]=un;break;
  1914.     case 19: y=cgetg(3,17);y[2]=lstoi(lg(x)-1);
  1915.       y[1]=(lg(x)>1)?lstoi(lg(x[1])-1):zero;break;
  1916.     default: err(matler1);
  1917.     }
  1918.   return y;
  1919. }
  1920.  
  1921. GEN geval(x)
  1922.   GEN x;
  1923. {
  1924.   long av, tetpil, tx = typ(x), lx, i;
  1925.   GEN y, z;
  1926.   if (tx < 9) return gcopy(x);
  1927.   if (tx > 13)
  1928.   {
  1929.     lx = lg(x);
  1930.     y = cgetg(lx, tx);
  1931.     for(i = 1; i < lx; i++) y[i] = (long)geval(x[i]);
  1932.     return y;
  1933.   }
  1934.   switch(tx)
  1935.   {
  1936.     case 9:
  1937.       y=cgetg(3,9);y[1]=(long)geval(x[1]);av=avma;
  1938.       z=geval(x[2]);tetpil=avma;y[2]=lpile(av,tetpil,gmod(z,y[1]));
  1939.       return y;
  1940.     case 10:
  1941.       lx = lgef(x); if (lx == 2) return gzero;
  1942.       y = gzero; z = (GEN)varentries[varn(x)]->value;
  1943.       av = avma;
  1944.       for(i = lx-1; i > 1; i--)
  1945.       {
  1946.         tetpil = avma;
  1947.         y = gadd(geval(x[i]), gmul(z, y));
  1948.       }
  1949.       return gerepile(av, tetpil, y);
  1950.     case 11: err(impl, "evaluation of a power series");
  1951.     case 13: return gdiv(geval(x[1]),geval(x[2]));
  1952.   }
  1953. }
  1954.  
  1955. int isexactzero(g)
  1956.      GEN g;
  1957. {
  1958.   long i;
  1959.   switch (typ(g))
  1960.   {
  1961.     case 1: return !signe(g);
  1962.     case 2:
  1963.     case 7:
  1964.     case 11: return 0;
  1965.     case 3:
  1966.     case 9: return isexactzero(g[2]);
  1967.     case 4:
  1968.     case 5:
  1969.     case 13:
  1970.     case 14: return isexactzero(g[1]);
  1971.     case 6: return isexactzero(g[1])&&isexactzero(g[2]);
  1972.     case 8: return isexactzero(g[2])&&isexactzero(g[3]);
  1973.     case 10: for (i=lgef(g)-1;i>1;i--) if (!isexactzero(g[i])) return 0;
  1974.       return 1;
  1975.     case 17:
  1976.     case 18:
  1977.     case 19: for(i=1;i<lg(g);i++) if(!isexactzero(g[i])) return 0;
  1978.       return 1;
  1979.     default: return 0;
  1980.   }
  1981. }
  1982.  
  1983. GEN simplify(x)
  1984.      GEN x;
  1985. {
  1986.   long tx=typ(x),av,tetpil,i,lx;
  1987.   GEN p1,y;
  1988.  
  1989.   switch(tx)
  1990.     {
  1991.     case 1:
  1992.     case 2:
  1993.     case 3:
  1994.     case 4:
  1995.     case 7:
  1996.     case 15:
  1997.     case 16: return gcopy(x);
  1998.     case 5: return gred(x);
  1999.     case 6:
  2000.     case 8: if(isexactzero(x[2])) return simplify(x[1]);
  2001.       else 
  2002.     {
  2003.       y=cgetg(lg(x),tx);y[1]=(long)simplify(x[1]);
  2004.       y[2]=(long)simplify(x[2]);return y;
  2005.     }
  2006.     case 9: y=cgetg(3,9);y[1]=(long)simplify(x[1]);
  2007.       av=avma;p1=gmod(x[2],y[1]);tetpil=avma;
  2008.       y[2]=lpile(av,tetpil,simplify(p1));return y;
  2009.     case 10: lx=lgef(x);if(lx==2) return gzero;
  2010.       if(lx==3) return simplify(x[2]);
  2011.       y=cgetg(lx,tx);y[1]=x[1];for(i=2;i<lx;i++) y[i]=(long)simplify(x[i]);
  2012.       return y;
  2013.     case 11: if (!signe(x)) return gcopy(x);
  2014.       lx=lg(x);
  2015.       y=cgetg(lx,tx);y[1]=x[1];for(i=2;i<lx;i++) y[i]=(long)simplify(x[i]);
  2016.       return y;
  2017.     case 13: y=cgetg(3,13);y[1]=(long)simplify(x[1]);
  2018.       y[2]=(long)simplify(x[2]);return y;
  2019.     case 14: av=avma;y=gred(x);tetpil=avma;
  2020.       return gerepile(av,tetpil,simplify(y));
  2021.     case 17:
  2022.     case 18:
  2023.     case 19: lx=lg(x);y=cgetg(lx,tx);
  2024.       for(i=1;i<lx;i++) y[i]=(long)simplify(x[i]);
  2025.       return y;
  2026.     default: err(simplifyer1);
  2027.     }
  2028. }
  2029.  
  2030. /* Karatsuba pour les polynomes */
  2031.  
  2032. GEN gmulxn(p,n,r)
  2033.      GEN p,*r;
  2034.      long n;
  2035.  
  2036. /* interne: pas de verifications et objects non connexes */
  2037.  
  2038. {
  2039.   long i,j,lx;
  2040.   GEN y;
  2041.  
  2042.   if((!signe(p))||(!n)) {*r=gzero;return gcopy(p);}
  2043.   lx=lgef(p);
  2044.   if(n>0)
  2045.     {
  2046.       y=cgetg(lx+n,10);y[1]=p[1]+n;
  2047.       for(i=2;i<=n+1;i++) y[i]=zero;
  2048.       for(i=n+2;i<lx+n;i++) y[i]=p[i-n];
  2049.       return y;
  2050.     }
  2051.   else
  2052.     {
  2053.       n= -n;if(n>lx-3) {*r=p;return gzero;}
  2054.       else
  2055.     {
  2056.       i=n+1;while((i>=2)&&gcmp0(p[i])) i--;
  2057.       if(i==1) *r=gzero;
  2058.       else 
  2059.         {
  2060.           y=cgetg(i+1,10);y[1]=p[1];setlgef(y,i+1);
  2061.           for(j=2;j<=i;j++) y[j]=p[j];*r=y;
  2062.         }
  2063.       y=cgetg(lx-n,10);y[1]=p[1]-n;
  2064.       for(i=2;i<lx-n;i++) y[i]=p[i+n];
  2065.       return y;
  2066.     }
  2067.     }
  2068. }
  2069.  
  2070. GEN karamul(x,y,k)
  2071.      GEN x,y;
  2072.      long k;
  2073. {
  2074.   long av=avma,tetpil,nx,ny,n;
  2075.   GEN p1,p2,p3,p4,p5,p6,p7;
  2076.  
  2077.   if((typ(x)!=10)||(typ(y)!=10)) err(karaer1);
  2078.   if(varn(x)!=varn(y)) return gmul(x,y);
  2079.   if((!signe(x))||(!signe(y))) {p1=cgetg(2,10);p1[1]=2;return p1;}
  2080.   if(k<=0) return gmul(x,y);
  2081.   nx=lgef(x)-3;ny=lgef(y)-3;n=(max(nx,ny)+1)>>1;
  2082.   p1=gmulxn(x,-n,&p2);p3=gmulxn(y,-n,&p4);
  2083.   p5=karamul(p1,p3,k-1);p6=karamul(p2,p4,k-1);
  2084.   p7=karamul(gadd(p1,p2),gadd(p3,p4),k-1);
  2085.   p1=gadd(gmulxn(p5,n<<1),gmulxn(gsub(p7,gadd(p5,p6)),n));
  2086.   tetpil=avma;return gerepile(av,tetpil,gadd(p1,p6));
  2087. }
  2088.   
  2089. /* karatsuba entier */
  2090.  
  2091. GEN mulxn(x,n,r)
  2092.      GEN x,*r;
  2093.      long n;
  2094. {
  2095.   long i,j,lx;
  2096.   GEN p1,p2;
  2097.  
  2098.   if((!n)||(!signe(x))) {*r=gzero;return x;}
  2099.   lx=lgef(x);
  2100.   if(n>0)
  2101.     {
  2102.       *r=gzero;p1=cgeti(lx+n);
  2103.       p1[1]=x[1]+n;for(i=2;i<lx;i++) p1[i]=x[i];
  2104.       for(i=lx;i<lx+n;i++) p1[i]=0;
  2105.       return p1;
  2106.     }
  2107.   else
  2108.     {
  2109.       n= -n;
  2110.       if(lx<=n+2) {*r=x;return gzero;}
  2111.       else
  2112.     {
  2113.       p1=cgeti(lx-n);p1[1]=x[1]-n;
  2114.       for(i=2;i<lx-n;i++) p1[i]=x[i];
  2115.       while((i<lx)&&(!x[i])) i++;
  2116.       if(i==lx) *r=gzero;
  2117.       else 
  2118.         {
  2119.           p2=cgeti(lx-i+2);p2[1]=0x01000002+lx-i;
  2120.           for(j=2;j<=lx-i+1;j++) p2[j]=x[j+i-2];
  2121.           *r=p2;
  2122.         }
  2123.       return p1;
  2124.     }
  2125.     }
  2126. }
  2127.  
  2128. GEN mpkaramul(x,y,k)
  2129.      GEN x,y;
  2130.      long k;
  2131. {
  2132.   long av=avma,tetpil,lx,ly,n;
  2133.   GEN x1,x2,y1,y2,p1,p2,p3;
  2134.  
  2135.   if((typ(x)!=1)||(typ(y)!=1)) err(karaer2);
  2136.   lx=lgef(x)-2;ly=lgef(y)-2;
  2137.   if((lx<=2)||(ly<=2)||(k<=0)) return mulii(x,y);
  2138.   n=(max(lx,ly)+1)>>1;
  2139.   x1=mulxn(x,-n,&x2);y1=mulxn(y,-n,&y2);
  2140.   p1=mpkaramul(x1,y1,k-1);p2=mpkaramul(x2,y2,k-1);
  2141.   p3=subii(mpkaramul(addii(x1,x2),addii(y1,y2),k-1),addii(p1,p2));
  2142.   p1=addii(mulxn(p1,n<<1,&x2),mulxn(p3,n,&x2));tetpil=avma;
  2143.   return gerepile(av,tetpil,addii(p1,p2));
  2144. }
  2145.  
  2146.   
  2147.   
  2148.