home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 2 / RISC_DISC_2.iso / pd_share / program / code / pari / c / gen3 < prev    next >
Encoding:
Text File  |  1991-11-28  |  48.2 KB  |  2,022 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. long    gvar(x)
  75.   
  76.      GEN  x;
  77.   
  78. {
  79.   long tx=typ(x),i,v,w;
  80.   if((tx<9)||(tx==15)||(tx==16)) return BIGINT;
  81.   if(tx==9) return varn(x[1]);
  82.   if(tx<12) return varn(x);
  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. long 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. 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. long 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.         if(ty==6) z=poleval(x,y);
  810.         else
  811.         {
  812.           vx=varn(x);l=lgef(x);
  813.           if(ty!=19)
  814.           {
  815.             if(vx>v) z=gcopy(x);
  816.             else
  817.             {
  818.               if(vx==v)
  819.               {
  820.                 if(l==3) z=gcopy(x[2]);
  821.                 else
  822.                 {
  823.                   av=avma;z=(GEN)(x[l-1]);
  824.                   for (i=l-1;i>3;i--)
  825.                     z=gadd(gmul(z,y),x[i-1]);
  826.                   z=gmul(z,y);tetpil=avma;
  827.                   z=gerepile(av,tetpil,gadd(z,x[2]));
  828.                 }
  829.               }
  830.               else
  831.               {
  832.                 if(l==3) z=gsubst(x[2],v,y);
  833.                 else
  834.                 {
  835.                   av=avma;p1=polx[vx];
  836.                   z= gsubst(x[l-1],v,y);
  837.                   for (i=l-1;i>3;i--)
  838.                     z=gadd(gmul(z,p1),gsubst(x[i-1],v,y));
  839.                   z=gmul(z,p1);p2=gsubst(x[2],v,y);
  840.                   tetpil=avma;
  841.                   z=gerepile(av,tetpil,gadd(z,p2));
  842.                 }
  843.               }
  844.             }
  845.           }
  846.           else
  847.           {
  848.             if(ly!=lg(y[1])) err(gsubser3);
  849.             if(vx>v) z=gscalmat(x,ly-1);
  850.             else
  851.             {
  852.               if(vx==v)
  853.               {
  854.                 if(l==3) z=gscalmat(x[2],ly-1);
  855.                 else
  856.                 {
  857.                   av=avma;z=(GEN)(x[l-1]);
  858.                   for (i=l-1;i>3;i--)
  859.                     z=gaddmat(x[i-1],gmul(z,y));
  860.                   z=gmul(z,y);tetpil=avma;
  861.                   z=gerepile(av,tetpil,gaddmat(x[2],z));
  862.                 }
  863.               }
  864.               else
  865.               {
  866.                 if(l==3) z=gsubst(x[2],v,y);
  867.                 else
  868.                 {
  869.                   av=avma;p1=polx[vx];
  870.                   z=gsubst(x[l-1],v,y);
  871.                   for(i=l-1;i>3;i--)
  872.                     z=gadd(gmul(z,p1),gsubst(x[i-1],v,y));
  873.                   z=gmul(z,p1);p2=gsubst(x[2],v,y);
  874.                   tetpil=avma;
  875.                   z=gerepile(av,tetpil,gadd(p2,z));
  876.                 }
  877.               }
  878.             }
  879.           }
  880.         }
  881.       }
  882.       break;
  883.       
  884.     case 11:
  885.       ex=valp(x);vx=varn(x);
  886.       if(vx>v)
  887.       {
  888.         if(ty==19) z=gscalmat(x,ly-1);
  889.         else z=gcopy(x);
  890.         return z;
  891.       }
  892.       if(!signe(x))
  893.       {
  894.         if(vx<v) z=gcopy(x);
  895.         else
  896.         {
  897.           z=cgetg(3,11);
  898.           z[1]=0x8000+ex*gval(y,v);z[2]=zero;
  899.           setvarn(z,vx);
  900.         }
  901.         return z;
  902.       }
  903.       if(vx<v)
  904.       {
  905.         /* a ameliorer */
  906.         av1=avma;setvalp(x,0);p1=gconvsp(x);setvalp(x,ex);
  907.         p2=gsubst(p1,v,y);av2=avma;
  908.         z=tayl(p2,vx,lx-2);
  909.         if(ex)
  910.         {
  911.           p1=gpuigs(polx[vx],ex);
  912.           av2=avma;z=gerepile(av1,av2,gmul(z,p1));
  913.         }
  914.         else z=gerepile(av1,av2,z);
  915.         return z;
  916.       }
  917.       switch(ty)
  918.       {
  919.         case 11:
  920.     ey=valp(y);vy=varn(y);
  921.     if (ey<1)
  922.       {
  923.         z=cgetg(3,11);z[2]=zero;
  924.         z[1]=0x8000+ey*(ex+lx-2);
  925.         setvarn(z,vy);
  926.       }
  927.     else
  928.       {
  929.         l=(lx-2)*ey+2;
  930.         if (ex)
  931.           {if (l>ly) l=ly;}
  932.         else 
  933.           {
  934.         if (gcmp0(y)) l=ey+2;
  935.         else
  936.           {if (l>ly) l=ly+ey;}
  937.           }
  938.         if(vy!=vx)
  939.           {
  940.         av=avma;z=cgetg(2,11);z[1]=0x8000;setvarn(z,vy);
  941.         for(i=lx-1;i>=2;i--)
  942.           {p1=gmul(y,z);av2=avma;z=gadd(x[i],p1);}
  943.         if (ex)
  944.           {
  945.             p1=gpuigs(y,ex);av2=avma;
  946.             z=gerepile(av,av2,gmul(z,p1));
  947.           }
  948.         else z=gerepile(av,av2,z);
  949.           }
  950.         else
  951.           {
  952.         av1=avma;t=cgetg(ly,11);
  953.         av2=avma;z=cgetg(l,11);
  954.         z[2] =lcopy(x[2]);
  955.         for (i=3;i<l;i++) z[i]=zero;
  956.         for (i=2;i<ly;i++) t[i]=y[i];
  957.         
  958.         for (i=3,jb=ey;jb<=l-2;i++,jb+=ey)
  959.           {
  960.             for (j=jb+2;j<l;j++)
  961.               {
  962.             av4=avma;
  963.             p1=gmul(x[i],t[j-jb]);
  964.             av5=avma;
  965.             z[j]=lpile(av4,av5,ladd(z[j],p1));
  966.             if (j==jb+ey+1) av=avma;
  967.               }
  968.             for (j=l-1-jb-ey;j>1;j--)
  969.               {
  970.             av4=avma;p1=gzero;
  971.             for (k=2;k<j;k++)
  972.               p1=gadd(p1,gmul(t[j-k+2],y[k]));
  973.             p2=gmul(t[2],y[j]);
  974.             av5=avma;
  975.             t[j]=lpile(av4,av5,gadd(p1,p2));
  976.               }
  977.             if (i>3)
  978.               {
  979.             av7=avma;decal=lpile(av3,av6,0);
  980.             for (k=jb+2;k<l;k++)
  981.               if (adecaler(z[k],av6,av7)) z[k]+=decal;
  982.             for (k=2;k<l-jb-ey+2;k++)
  983.               if (adecaler(t[k],av6,av7)) t[k]+=decal;
  984.             av3=av+decal;
  985.               }
  986.             else av3=av;
  987.             av6=avma;
  988.           }
  989.         z[1]=0x1008000;setvarn(z,varn(y));
  990.         if (ex)
  991.           {
  992.             if (l<ly) setlg(y,l);
  993.             p1=gpuigs(y,ex);
  994.             av2=avma;
  995.             z=gerepile(av1,av2,gmul(z,p1));
  996.             if (l<ly) setlg(y,ly);
  997.           }
  998.         else z=gerepile(av1,av2,z);
  999.           }
  1000.       }
  1001.           break;
  1002.           
  1003.         case 10:
  1004.         case 13:
  1005.         case 14:
  1006.           if(gcmp0(y))
  1007.           {
  1008.             z=cgetg(lx,tx);z[1]=0x1008000;
  1009.             z[2]=lcopy(x[2]);setvarn(z,v);
  1010.             for(i=3;i<lx;i++) z[i]=zero;
  1011.           }
  1012.           else
  1013.           {
  1014.             e=gval(y,vy=gvar(y));if(e<=0) err(gsubser5);
  1015.         av=avma;p1=gconvsp(x);p2=gsubst(p1,v,y);tetpil=avma;
  1016.         z=gerepile(av,tetpil,tayl(p2,vy,e*(lx-2)));
  1017.           }
  1018.           break;
  1019.         default: err(gsubser4);
  1020.         }
  1021.       break;
  1022.     case 9:
  1023.       z=cgetg(lx,tx);z[1]=lsubst(x[1],v,y);av=avma;
  1024.       p1=gsubst(x[2],v,y);tetpil=avma;
  1025.       z[2]=lpile(av,tetpil,gmod(p1,z[1]));
  1026.       break;
  1027. /*  case 9: err(gsubser1); */
  1028.     case 13:
  1029.     case 14:
  1030.       av=avma;p1=gsubst(x[1],v,y);p2=gsubst(x[2],v,y);
  1031.       tetpil=avma;z=gerepile(av,tetpil,gdiv(p1,p2));
  1032.       break;
  1033.     case 17:
  1034.     case 18:
  1035.     case 19:
  1036.       z=cgetg(lx,tx);
  1037.       for(i=1;i<lx;i++)
  1038.       z[i]=lsubst(x[i],v,y);
  1039.   }
  1040.   return z;
  1041. }
  1042.  
  1043. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1044. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1045. /*                                                                 */
  1046. /*                SERIE RECIPROQUE D'UNE SERIE                     */
  1047. /*                                                                 */
  1048. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1049. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1050.  
  1051. GEN     recip(x)
  1052.      
  1053.      GEN     x;
  1054.      
  1055. {
  1056.   
  1057.   long  lx,av1,av2,av3,av4,av5,av6,av7,av8,i,j,k,v,decal;
  1058.   GEN   a,y,u,p1,p2;
  1059.   
  1060.   if ((typ(x)!=11) || (valp(x)!=1)) err(reciper);
  1061.   /* attention au coeff directeur */
  1062.   a=(GEN)x[2];v=varn(x);
  1063.   if (gcmp1(a))
  1064.   {
  1065.     lx=lg(x);av1=avma;u=cgetg(lx,11);
  1066.     av2=avma;y=cgetg(lx,11);
  1067.     u[2]=un;y[2]=un;
  1068.     av5=avma;u[3]=lmulsg(-2,x[3]);
  1069.     av6=avma;y[3]=lneg(x[3]);
  1070.     av7=avma;i=2;
  1071.     while (++i<lx-1)
  1072.     {
  1073.       for (j=3;j<i+1;j++)
  1074.       {
  1075.         av3=avma;p1=(GEN)u[j];
  1076.         for (k=j-1;k>2;k--)
  1077.           p1=gsub(p1,gmul(u[k],x[j-k+2]));
  1078.         av4=avma;
  1079.         u[j]=lpile(av3,av4,lsub(p1,x[j]));
  1080.       }
  1081.       av3=avma;p1=gmulsg(i,x[i+1]);
  1082.       for (k=2;k<i;k++)
  1083.       {
  1084.         p2=gmul(x[k+1],u[i-k+2]);
  1085.         p1=gadd(p1,gmulsg(k,p2) );
  1086.       }
  1087.       av4=avma;u[i+1]=lpile(av3,av4,lneg(p1));
  1088.       av8=avma;decal=lpile(av5,av6,0);
  1089.       for (k=3;k<i+2;k++)
  1090.         if (adecaler(u[k],av6,av8)) u[k]+=decal;
  1091.       if(adecaler(y[i],av6,av8)) y[i]+=decal;
  1092.       av6=avma;av5=av7+decal;
  1093.       y[i+1]=ldivgs(u[i+1],i);av7=avma;
  1094.     }
  1095.     y[lx-1]=lpile(av5,av6,y[lx-1]);
  1096.     y[1]=0x1008001;y=gerepile(av1,av2,y);
  1097.     setvarn(y,v);
  1098.   }
  1099.   else
  1100.   {
  1101.     p1=(GEN)x[2];av1=avma;y=gdiv(x,p1);
  1102.     y=recip(y);p2=gdiv(polx[v],p1);av2=avma;
  1103.     y=gerepile(av1,av2,gsubst(y,v,p2));
  1104.   }
  1105.   return y;
  1106. }
  1107.  
  1108.  
  1109. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1110. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1111. /*                                                                 */
  1112. /*                    DERIVATION ET INTEGRATION                    */
  1113. /*                                                                 */
  1114. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1115. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1116.  
  1117. GEN     deriv(x,v)
  1118.      
  1119.      GEN     x;
  1120.      long    v;
  1121.      
  1122. {
  1123.   long  lx,ly,vx,tx,e,i,j,tetpil,l,l1,f;
  1124.   GEN   y,p1,p2;
  1125.   
  1126.   tx=typ(x);if(tx<9) return gzero;
  1127.   switch(tx)
  1128.     {
  1129.     case 9: if(v<=varn(x[1])) return gzero;
  1130.       y=cgetg(3,9);y[1]=copyifstack(x[1]);y[2]=lderiv(x[2],v);
  1131.       return y;break;
  1132.     case 10: vx=varn(x);
  1133.       if(vx>v) return gzero;
  1134.       lx=lgef(x)-1;
  1135.       if(vx<v)
  1136.       {
  1137.         l=avma;
  1138.         for(i=lx;(i>=2)&&(isexactzero(p1=deriv(x[i],v)));avma=l,i--);
  1139.         y=cgetg(i+1,10);
  1140.         if(i==1) y[1]=2;
  1141.         else
  1142.         {
  1143.           y[1]=0x1000001+i;y[i]=(long)p1;f=gcmp0(p1);
  1144.           for(j=2;j<i;j++) {y[j]=lderiv(x[j],v);if(f) f=gcmp0(y[j]);}
  1145.       if(f) setsigne(y,0);
  1146.         }
  1147.         setvarn(y,vx);
  1148.         return y;
  1149.       }
  1150.       if(lx<3) y=gzero;
  1151.       else
  1152.       {
  1153.         l=avma;
  1154.         for(i=lx-1;(i>=2)&&isexactzero(p1=gmulsg(i-1,x[i+1]));avma=l,i--);
  1155.         y=cgetg(i+1,tx);
  1156.         if(i==1) y[1]=2;
  1157.         else
  1158.         {
  1159.           y[1]=0x1000001+i;y[i]=(long)p1;f=gcmp0(p1);
  1160.           for(i--;i>=2;i--) {y[i]=lmulsg(i-1,x[i+1]);if(f) f=gcmp0(y[i]);}
  1161.       if(f) setsigne(y,0);
  1162.         }
  1163.         setvarn(y,v);
  1164.       }
  1165.       break;
  1166.     case 11: vx=varn(x);
  1167.       if(vx>v) return gzero;
  1168.       lx=lg(x);e=valp(x);
  1169.       if(vx<v)
  1170.       {
  1171.         if(!signe(x)) y=gcopy(x);
  1172.         else
  1173.         {
  1174.           l=avma;
  1175.           for(i=2;(i<lx)&&(gcmp0(p1=deriv(x[i],v)));avma=l,i++);
  1176.           if(i==lx) y=ggrando(polx[vx],e+lx-2);
  1177.           else
  1178.           {
  1179.             y=cgetg(lx-i+2,11);y[1]=0x7ffe +e+i;
  1180.             setvarn(y,vx);y[2]=(long)p1;
  1181.             for(j=3;j<=lx-i+1;j++)
  1182.               y[j]=lderiv(x[i+j-2],v);
  1183.           }
  1184.         }
  1185.         return y;
  1186.       }
  1187.       ly=lx-1;if(ly<3) ly=3;
  1188.       if(gcmp0(x))
  1189.         {y=cgetg(3,tx);y[1]=0x7fff+e;}
  1190.       else
  1191.       {
  1192.         if(e)
  1193.         {
  1194.           y=cgetg(lx,tx);
  1195.           for(i=2;i<lx;i++)
  1196.             y[i]=lmulsg(i+e-2,x[i]);
  1197.           y[1]=0x1008000+e-1;
  1198.         }
  1199.         else
  1200.         {
  1201.           i=3;
  1202.           while((i<lx)&&gcmp0(x[i])) i++;
  1203.           if(i==lx)
  1204.             {y=cgetg(ly,tx);y[1]=0x7ffd+lx;}
  1205.           else
  1206.           {
  1207.             ly=ly-i+3;y=cgetg(ly,tx);
  1208.             y[1]=0x1007ffd+i;
  1209.             for(j=2;j<ly;j++)
  1210.               y[j]=lmulsg(j+i-4,x[i+j-2]);
  1211.           }
  1212.         }
  1213.       }
  1214.       setvarn(y,vx);
  1215.       break;
  1216.     case 13:
  1217.     case 14: l=avma;y=cgetg(3,tx);y[2]=lmul(x[2],x[2]);
  1218.       l1=avma;p1=gmul(x[2],deriv(x[1],v));p2=gmul(x[1],deriv(x[2],v));
  1219.       tetpil=avma;p1=gsub(p1,p2);y[1]=lpile(l1,tetpil,p1);
  1220.       if(tx==13) {tetpil=avma;y=gerepile(l,tetpil,gred(y));}
  1221.       break;
  1222.     case 17:
  1223.     case 18:
  1224.     case 19: lx=lg(x);y=cgetg(lx,tx);
  1225.       for(i=1;i<lx;i++)
  1226.         y[i]=lderiv(x[i],v);
  1227.       break;
  1228.     default: break;
  1229.     }
  1230.   return y;
  1231. }
  1232.  
  1233.  
  1234. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1235. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1236. /*                                                                 */
  1237. /*                    INTEGRATION FORMELLE                         */
  1238. /*                                                                 */
  1239. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1240. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1241.  
  1242. GEN     integ(x,v)
  1243.      
  1244.      GEN     x;
  1245.      long    v;
  1246.      
  1247. {
  1248.   long  lx,tx,e,i,j,vx;
  1249.   GEN   y;
  1250.   
  1251.   tx=typ(x);
  1252.   if((tx<9)||((tx==9)&&(v<=varn(x[1]))))
  1253.   {
  1254.     if(gcmp0(x)) y=gzero;
  1255.     else
  1256.     {
  1257.       y=cgetg(4,10);y[1]=0x1000004;setvarn(y,v);
  1258.       y[2]=zero;y[3]=lcopy(x);
  1259.     }
  1260.     return y;
  1261.   }
  1262.   switch(tx)
  1263.     {
  1264.     case 9: y=cgetg(3,9);y[1]=copyifstack(x[1]);y[2]=linteg(x[2],v);
  1265.       return y;break;
  1266.     case 10: vx=varn(x);lx=lgef(x);
  1267.       if(lx==2)
  1268.       {
  1269.         y=cgetg(3,10);y[1]=2;
  1270.         if(vx>v) setvarn(y,v);else setvarn(y,vx);
  1271.         return y;
  1272.       }
  1273.       if(vx>v)
  1274.       {
  1275.         y=cgetg(4,10);y[1]=(signe(x))?0x1000004:4;setvarn(y,v);
  1276.         y[2]=zero;y[3]=lcopy(x);
  1277.         return y;
  1278.       }
  1279.       if(vx<v)
  1280.       {
  1281.         y=cgetg(lx,10);y[1]=x[1];
  1282.         for(i=2;i<lx;i++)
  1283.           y[i]=linteg(x[i],v);
  1284.         return y;
  1285.       }
  1286.       y=cgetg(lx+1,tx);y[2]=zero;
  1287.       for(i=3;i<=lx;i++)
  1288.         y[i]=ldivgs(x[i-1],i-2);
  1289.       y[1]=signe(x)?0x1000001+lx:1+lx;setvarn(y,v);
  1290.       break;
  1291.     case 11: lx=lg(x);e=valp(x);vx=varn(x);
  1292.       if(!signe(x))
  1293.       {
  1294.         y=cgetg(3,tx);
  1295.         if(vx==v) y[1]=0x8001+e;
  1296.         else y[1]=x[1];
  1297.         if(vx>v) setvarn(y,v);
  1298.         return y;
  1299.       }
  1300.       if(vx>v)
  1301.       {
  1302.         y=cgetg(4,10);y[1]=0x1000004;setvarn(y,v);
  1303.         y[2]=zero;y[3]=lcopy(x);
  1304.         return y;
  1305.       }
  1306.       if(vx<v)
  1307.       {
  1308.         y=cgetg(lx,tx);y[1]=x[1];
  1309.         for(i=2;i<lx;i++)
  1310.           y[i]=linteg(x[i],v);
  1311.         return y;
  1312.       }
  1313.       y=cgetg(lx,tx);
  1314.       for(i=2;i<lx;i++)
  1315.       {
  1316.         if(!(j=i+e-1)) err(inter2);
  1317.         y[i]=ldivgs(x[i],j);
  1318.       }
  1319.       y[1]=(x[1])+1;
  1320.       break;
  1321.     case 13:
  1322.     case 14: err(impl,"integration of a rational function");
  1323.     case 17:
  1324.     case 18:
  1325.     case 19: lx=lg(x);y=cgetg(lx,tx);
  1326.       for(i=1;i<lx;i++)
  1327.         y[i]=linteg(x[i],v);
  1328.       break;
  1329.     default: err(inter1);
  1330.   }
  1331.   return y;
  1332. }
  1333.  
  1334. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1335. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1336. /*                                                                 */
  1337. /*                    PARTIES ENTIERES                             */
  1338. /*                                                                 */
  1339. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1340. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1341.  
  1342. GEN     gfloor(x)
  1343.      
  1344.      GEN     x;
  1345.      
  1346. {
  1347.   GEN     y,p1;
  1348.   long  i,lx,tx,av,tetpil;
  1349.   
  1350.   switch(tx=typ(x))
  1351.   {
  1352.     case 1 :
  1353.     case 10: y=gcopy(x);break;
  1354.     case 2 : y=mpent(x);break;
  1355.     case 4 :
  1356.     case 5 : av=avma;y=dvmdii(x[1],x[2],&p1);
  1357.       i=!gcmp0(p1);cgiv(p1);
  1358.       if(i&&(gsigne(x)<0))
  1359.       {
  1360.         tetpil=avma;y=gerepile(av,tetpil,gsubgs(y,1));
  1361.       }
  1362.       break;
  1363.     case 13:
  1364.     case 14: y=gdeuc(x[1],x[2]);
  1365.       break;
  1366.     case 17:
  1367.     case 18:
  1368.     case 19: lx=lg(x);y=cgetg(lx,tx);
  1369.       for(i=1;i<lx;i++)
  1370.         y[i]=lfloor(x[i]);
  1371.       break;
  1372.     default: err(flooer);
  1373.   }
  1374.   return y;
  1375. }
  1376.  
  1377. GEN     gfrac(x)
  1378.      
  1379.      GEN        x;
  1380.      
  1381. {
  1382.   long  av,tetpil;
  1383.   GEN   p1;
  1384.   
  1385.   av=avma;p1=gfloor(x);tetpil=avma;
  1386.   return gerepile(av,tetpil,gsub(x,p1));
  1387. }
  1388.  
  1389. GEN     gceil(x)
  1390.      
  1391.      GEN     x;
  1392.      
  1393. {
  1394.   GEN     y,p1;
  1395.   long  i,lx,tx=typ(x),av,tetpil;
  1396.   
  1397.   switch(tx)
  1398.   {
  1399.     case 1 :
  1400.     case 10: y=gcopy(x);break;
  1401.     case 2 : av=avma;y=mpent(x);
  1402.       if(!gegal(x,y))
  1403.       {
  1404.         tetpil=avma;y=gerepile(av,tetpil,gaddsg(1,y));
  1405.       }
  1406.       break;
  1407.     case 4 :
  1408.     case 5 : av=avma;y=dvmdii(x[1],x[2],&p1);
  1409.       i=!gcmp0(p1);cgiv(p1);
  1410.       if(i&&(gsigne(x)>0))
  1411.       {
  1412.         tetpil=avma;y=gerepile(av,tetpil,gaddsg(1,y));
  1413.       }
  1414.       break;
  1415.     case 13:
  1416.     case 14: y=gdeuc(x[1],x[2]);
  1417.       break;
  1418.     case 17:
  1419.     case 18:
  1420.     case 19: lx=lg(x);y=cgetg(lx,tx);
  1421.       for(i=1;i<lx;i++)
  1422.         y[i]=lceil(x[i]);
  1423.       break;
  1424.     default: err(ceiler);
  1425.   }
  1426.   return y;
  1427. }
  1428.  
  1429. GEN     ground(x)
  1430.      
  1431.      GEN     x;
  1432.      
  1433. {
  1434.   GEN     y,p1;
  1435.   long  i,lx=lg(x),tx=typ(x),av,tetpil;
  1436.   
  1437.   switch(tx)
  1438.   {
  1439.     case 1 :
  1440.     case 3 :
  1441.     case 8 : y=gcopy(x);break;
  1442.     case 2 : 
  1443.     case 4 :
  1444.     case 5 : av=avma;p1=gadd(ghalf,x);tetpil=avma;
  1445.       y=gerepile(av,tetpil,gfloor(p1));break;
  1446.     case 9 : y=cgetg(3,9);y[1]=copyifstack(x[1]);y[2]=lround(x[2]);
  1447.       break;
  1448.     case 10: lx=lgef(x);
  1449.     case 6 :
  1450.     case 11:
  1451.     case 13:
  1452.     case 14: 
  1453.     case 17:
  1454.     case 18:
  1455.     case 19: lx=lg(x);y=cgetg(lx,tx);
  1456.       for(i=1;i<lontyp[tx];i++) y[i]=x[i];
  1457.       for(i=lontyp[tx];i<lx;i++) y[i]=lround(x[i]);
  1458.       break;
  1459.     default: err(rounder);
  1460.   }
  1461.   return y;
  1462. }
  1463.  
  1464. GEN     grndtoi(x,e)
  1465.      
  1466.      GEN     x;
  1467.      long   *e;
  1468.      
  1469. {
  1470.   GEN     y,p1;
  1471.   long  i,lx=lg(x),tx=typ(x),av,tetpil,ex,e1;
  1472.   
  1473.   *e= -1048576;
  1474.   switch(tx)
  1475.   {
  1476.     case 1 :
  1477.     case 3 :
  1478.     case 8 : y=gcopy(x);break;
  1479.     case 2 : av=avma;p1=gadd(ghalf,x);
  1480.       ex=expo(p1);*e=ex-(lx-2)*32;
  1481.       if(ex<0)
  1482.       {
  1483.     if(signe(p1)>=0) {avma=av;y=gzero;}
  1484.     else {tetpil=avma;y=gerepile(av,tetpil,gneg(un));}
  1485.       }
  1486.       else
  1487.       {
  1488.         settyp(p1,1);setlgef(p1,lx);
  1489.         if(signe(x)>=0)
  1490.           {tetpil=avma;y=gerepile(av,tetpil,shifti(p1,*e+1));}
  1491.         else
  1492.         {
  1493.           y=shifti(p1,*e+1);tetpil=avma;
  1494.           y=gerepile(av,tetpil,gaddgs(y,-1));
  1495.         }
  1496.       }
  1497.       break;
  1498.     case 4 :
  1499.     case 5 : av=avma;p1=gadd(ghalf,x);tetpil=avma;
  1500.       y=gerepile(av,tetpil,gfloor(p1));break;
  1501.     case 9 : y=cgetg(3,9);y[1]=copyifstack(x[1]);y[2]=lrndtoi(x[2],e);
  1502.       break;
  1503.     case 10: lx=lgef(x);
  1504.     case 6 :
  1505.     case 11:
  1506.     case 13:
  1507.     case 14: 
  1508.     case 17:
  1509.     case 18:
  1510.     case 19: lx=lg(x);y=cgetg(lx,tx);
  1511.       for(i=1;i<lontyp[tx];i++)
  1512.         y[i]=x[i];
  1513.       for(i=lontyp[tx];i<lx;i++)
  1514.       {
  1515.         y[i]=lrndtoi(x[i],&e1);
  1516.         if(e1>*e) *e=e1;
  1517.       }
  1518.       break;
  1519.     default: err(rndtoier);
  1520.   }
  1521.   return y;
  1522. }
  1523.  
  1524. long rounderror(x)
  1525.      GEN x;
  1526.      
  1527. {
  1528.   long e,av=avma;
  1529.  
  1530.   grndtoi(x,&e);avma=av;
  1531.   e*=L2SL10;return e;
  1532. }
  1533.  
  1534. GEN     gcvtoi(x,e)
  1535.      
  1536.      GEN     x;
  1537.      long    *e;
  1538.      
  1539.      /* la variable formelle e,representant le nombre de bits d'erreur */
  1540.      /* sur la partie entiere,n'est utilisee que dans le cas 2 (reel) */
  1541.      
  1542. {
  1543.   GEN     y,p1;
  1544.   long  tx=typ(x),lx=lg(x),i,ex,av,tetpil,e1;
  1545.   
  1546.   *e= -1048576;
  1547.   switch(tx)
  1548.   {
  1549.     case 1 :
  1550.     case 10: y=gcopy(x);break;
  1551.     case 2 : ex=expo(x);*e=ex-(lx-2)*32;
  1552.       if(ex<0) y=gzero;
  1553.       else
  1554.       {
  1555.         av=avma;p1=gcopy(x);settyp(p1,1);setlgef(p1,lx);
  1556.         tetpil=avma;y=gerepile(av,tetpil,shifti(p1,*e+1));
  1557.       }
  1558.       break;
  1559.     case 4 :
  1560.     case 5 : y=divii(x[1],x[2]);  break;
  1561.     case 7 : y=gconvpe(x);break;
  1562.     case 13:
  1563.     case 14: y=gdeuc(x[1],x[2]);  break;
  1564.     case 11: y=gconvsp(x);break;
  1565.     case 17:
  1566.     case 18:
  1567.     case 19: y=cgetg(lx,tx);
  1568.       for(i=1;i<lx;i++)
  1569.       {
  1570.         y[i]=lcvtoi(x[i],&e1);
  1571.         if(e1>*e) *e=e1;
  1572.       }
  1573.       break;
  1574.     default: err(cvtoier);
  1575.   }
  1576.   return y;
  1577. }
  1578.  
  1579. GEN     gtrunc(x)
  1580.      
  1581.      GEN     x;
  1582.      
  1583. {
  1584.   GEN     y;
  1585.   long  tx=typ(x),lx=lg(x),i;
  1586.   switch(tx)
  1587.   {
  1588.     case 1 :
  1589.     case 10: y=gcopy(x);break;
  1590.     case 2 : y=mptrunc(x);break;
  1591.     case 4 :
  1592.     case 5 : y=divii(x[1],x[2]); break;
  1593.     case 7 : y=gconvpe(x);break;
  1594.     case 13:
  1595.     case 14: y=gdeuc(x[1],x[2]); break;
  1596.     case 11: y=gconvsp(x);break;
  1597.     case 17:
  1598.     case 18:
  1599.     case 19: y=cgetg(lx,tx);
  1600.       for(i=1;i<lx;i++)
  1601.         y[i]=ltrunc(x[i]);
  1602.       break;
  1603.     default: err(truncer);
  1604.   }
  1605.   return y;
  1606. }
  1607.  
  1608. GEN gtopoly(x,v)
  1609.      GEN x;
  1610.      long v;
  1611. {
  1612.   long tx=typ(x),lx,i,j;
  1613.   GEN y;
  1614.   
  1615.   if(gcmp0(x)) {y=cgetg(2,10);y[1]=2;setvarn(y,v);return y;}
  1616.   if(tx<10)
  1617.     {
  1618.       y=cgetg(3,10);y[1]=0x1000003;y[2]=lcopy(x);
  1619.       setvarn(y,v);return y;
  1620.     }
  1621.   switch(tx)
  1622.     {
  1623.     case 10: y=gcopy(x);break;
  1624.     case 11: y=gconvsp(x);break;
  1625.     case 13:
  1626.     case 14: y=gdeuc(x[1],x[2]);break;
  1627.     case 15:
  1628.     case 16:
  1629.     case 17:
  1630.     case 18:
  1631.     case 19: lx=lg(x);i=1;while((i<lx)&&gcmp0(x[i])) i++;
  1632.       y=cgetg(lx-i+2,10);y[1]=0x1000002+lx-i;
  1633.       for(j=2;j<=lx-i+1;j++) y[j]=lcopy(x[lx+1-j]);
  1634.     }
  1635.   setvarn(y,v);return y;
  1636. }
  1637.  
  1638. GEN gtopolyrev(x,v)
  1639.      GEN x;
  1640.      long v;
  1641. {
  1642.   long tx=typ(x),lx,i,j;
  1643.   GEN y;
  1644.   
  1645.   if(gcmp0(x)) {y=cgetg(2,10);y[1]=2;setvarn(y,v);return y;}
  1646.   if(tx<10)
  1647.     {
  1648.       y=cgetg(3,10);y[1]=0x1000003;y[2]=lcopy(x);
  1649.       setvarn(y,v);return y;
  1650.     }
  1651.   switch(tx)
  1652.     {
  1653.     case 10: y=gcopy(x);break;
  1654.     case 11: y=gconvsp(x);break;
  1655.     case 13:
  1656.     case 14: y=gdeuc(x[1],x[2]);break;
  1657.     case 15:
  1658.     case 16:
  1659.     case 17:
  1660.     case 18:
  1661.     case 19: lx=lg(x);i=1;while((i<lx)&&gcmp0(x[lx-i])) i++;
  1662.       y=cgetg(lx-i+2,10);y[1]=0x1000002+lx-i;
  1663.       for(j=2;j<=lx-i+1;j++) y[j]=lcopy(x[j-1]);
  1664.     }
  1665.   setvarn(y,v);return y;
  1666. }
  1667.  
  1668. GEN gtoser(x,v)
  1669.      GEN x;
  1670.      long v;
  1671. {
  1672.   long tx=typ(x),lx,i,j,l,av,tetpil;
  1673.   GEN y,p1,p2;
  1674.   
  1675.   if(tx==11) {y=gcopy(x);setvarn(y,v);return y;}
  1676.   if(gcmp0(x)) {y=cgetg(2,11);y[1]=0x8000+precdl;setvarn(y,v);return y;}
  1677.   if(tx<10)
  1678.     {
  1679.       y=cgetg(precdl+2,11);y[1]=0x1008000;y[2]=lcopy(x);
  1680.       for(i=3;i<=precdl+1;i++) y[i]=zero;
  1681.       setvarn(y,v);return y;
  1682.     }
  1683.   switch(tx)
  1684.     {
  1685.     case 10: lx=lgef(x);i=2;while((i<lx)&&gcmp0(x[i])) i++;
  1686.       l=lx-i;if(precdl>l) l=precdl;
  1687.       y=cgetg(l+2,11);y[1]=0x1008000+i-2;
  1688.       for(j=2;j<=lx-i+1;j++) y[j]=lcopy(x[j+i-2]);
  1689.       for(j=lx-i+2;j<=l+1;j++) y[j]=zero;
  1690.       break;
  1691.     case 13:
  1692.     case 14: av=avma;p1=gtoser(x[1],v);p2=gtoser(x[2],v);
  1693.       tetpil=avma;y=gerepile(av,tetpil,gdiv(p1,p2));break;
  1694.     case 15:
  1695.     case 16:
  1696.     case 17:
  1697.     case 18:
  1698.     case 19: lx=lg(x);i=1;while((i<lx)&&gcmp0(x[i])) i++;
  1699.       y=cgetg(lx-i+2,11);y[1]=0x1008000+i-1;
  1700.       for(j=2;j<=lx-i+1;j++) y[j]=lcopy(x[j+i-2]);
  1701.     }
  1702.   setvarn(y,v);return y;
  1703. }
  1704.  
  1705. GEN gtovec(x)
  1706.      GEN x;
  1707. {
  1708.   long tx=typ(x),lx,i;
  1709.   GEN y;
  1710.   
  1711.   if((tx<10)||(tx==13)||(tx==14))
  1712.     {y=cgetg(2,17);y[1]=lcopy(x);return y;}
  1713.   if(tx>=15)
  1714.     {
  1715.       lx=lg(x);y=cgetg(lx,17);
  1716.       for(i=1;i<lx;i++) y[i]=lcopy(x[i]);
  1717.       return y;
  1718.     }
  1719.   if(tx==10)
  1720.     {
  1721.       lx=lgef(x);y=cgetg(lx-1,17);
  1722.       for(i=1;i<=lx-2;i++) y[i]=lcopy(x[lx-i]);
  1723.       return y;
  1724.     }
  1725.   if(!signe(x)) {y=cgetg(2,17);y[1]=zero;return y;}
  1726.   lx=lg(x);y=cgetg(lx-1,17);
  1727.   for(i=1;i<=lx-2;i++) y[i]=lcopy(x[i+1]);
  1728.   return y;
  1729. }
  1730.       
  1731. GEN     compo(x,n)
  1732.      
  1733.      GEN     x;
  1734.      long    n;
  1735.      
  1736. {
  1737.   long l,lx=lg(x),tx=typ(x);
  1738.  
  1739.   if((tx==10)&&((n+1)>=lgef(x))) return gzero;
  1740.   if((tx==11)&&(!signe(x))) return gzero;
  1741.   l=lontyp[tx]+n-1;
  1742.   if((l>=lx) || (n<1)) err(compoer1);
  1743.   return gcopy(x[l]);
  1744. }
  1745.  
  1746. GEN truecoeff(x,n)
  1747.      GEN x;
  1748.      long n;
  1749. {
  1750.   long tx=typ(x),lx,ex;
  1751.  
  1752.   if(tx<10)
  1753.     {if(n) return gzero;else return gcopy(x);}
  1754.   switch(tx)
  1755.     {
  1756.     case 15: case 16: case 17: case 18: case 19:
  1757.       if((n<=0)||(n>=lg(x))) err(compoer1);
  1758.       return gcopy(x[n]);break;
  1759.     case 10:
  1760.       if((n<0)||(n>=lgef(x)-2)) return gzero;
  1761.       return gcopy(x[n+2]);break;
  1762.     case 11:
  1763.       if(!signe(x)) return gzero;
  1764.       lx=lg(x);ex=valp(x);
  1765.       if(n<ex) return gzero;
  1766.       if(n>=ex+lx-2) err(compoer1);
  1767.       return gcopy(x[n-ex+2]);break;
  1768.     default: err(compoer1);
  1769.     }
  1770. }
  1771.  
  1772. GEN  denom(x)
  1773.      GEN x;
  1774. {
  1775.   long i,av,tetpil;
  1776.   GEN s;
  1777.  
  1778.   switch(typ(x))
  1779.     {
  1780.     case 1: return gun;
  1781.     case 4: case 5: return absi(x[2]);
  1782.     case 13: case 14: return gcopy(x[2]);
  1783.     case 10: return polun[varn(x)];
  1784.     case 17: case 18: case 19: av=tetpil=avma;
  1785.       s=denom(x[1]);
  1786.       for(i=2;(i<lg(x));i++) {tetpil=avma;s=glcm(s,denom(x[i]));}
  1787.       return gerepile(av,tetpil,s);
  1788.     default: err(denomer1);
  1789.     }
  1790. }
  1791.  
  1792. GEN  numer(x)
  1793.      GEN x;
  1794. {
  1795.   long av,tetpil;
  1796.   GEN s;
  1797.  
  1798.   switch(typ(x))
  1799.     {
  1800.     case 1: case 10: return gcopy(x);
  1801.     case 4: case 5: return (signe(x[2])>0) ? gcopy(x[1]) : gneg(x[1]);
  1802.     case 13: case 14: return gcopy(x[1]);
  1803.     case 17: case 18: case 19: av=avma;s=denom(x);tetpil=avma;
  1804.       return gerepile(av,tetpil,gmul(s,x));
  1805.     default: err(numerer1);
  1806.     }
  1807. }
  1808.  
  1809. GEN lift(x)
  1810.      GEN x;
  1811. {
  1812.   long lx,tx=typ(x),i;
  1813.   GEN y;
  1814.  
  1815.   switch(tx)
  1816.     {
  1817.     case 1: return gcopy(x);
  1818.     case 3:
  1819.     case 9: return gcopy(x[2]);
  1820.     case 11: if(!signe(x)) return gcopy(x);
  1821.       y=cgetg(lx=lg(x),tx);y[1]=x[1];
  1822.       for(i=2;i<lx;i++) y[i]=(long)lift(x[i]);break;
  1823.     case 4:
  1824.     case 5:
  1825.     case 6:
  1826.     case 13:
  1827.     case 14:
  1828.     case 17:
  1829.     case 18:
  1830.     case 19: y=cgetg(lx=lg(x),tx);
  1831.       for(i=1;i<lx;i++) y[i]=(long)lift(x[i]);
  1832.       break;
  1833.     case 10: y=cgetg(lx=lgef(x),tx);y[1]=x[1];
  1834.       for(i=2;i<lx;i++) y[i]=(long)lift(x[i]);break;
  1835.     case 8: y=cgetg(4,tx);y[1]=copyifstack(x[1]);
  1836.       for(i=2;i<4;i++) y[i]=(long)lift(x[i]);break;
  1837.     default: err(lifter1);
  1838.     }
  1839.   return y;
  1840. }
  1841.  
  1842. GEN centerlift(x)
  1843.      GEN x;
  1844. {
  1845.   long lx,tx=typ(x),i,av;
  1846.   GEN y;
  1847.  
  1848.   switch(tx)
  1849.     {
  1850.     case 1: return gcopy(x);
  1851.     case 3: av=avma;i=cmpii(shifti(x[2],1),x[1]);avma=av;
  1852.       y=(i>0)?subii(x[2],x[1]):gcopy(x[2]);break;
  1853.     case 9: y=gcopy(x[2]);break;
  1854.     case 11: if(!signe(x)) return gcopy(x);
  1855.       y=cgetg(lx=lg(x),tx);y[1]=x[1];
  1856.       for(i=2;i<lx;i++) y[i]=(long)centerlift(x[i]);break;
  1857.     case 4:
  1858.     case 5:
  1859.     case 6:
  1860.     case 13:
  1861.     case 14:
  1862.     case 17:
  1863.     case 18:
  1864.     case 19: y=cgetg(lx=lg(x),tx);
  1865.       for(i=1;i<lx;i++) y[i]=(long)centerlift(x[i]);
  1866.       break;
  1867.     case 10: y=cgetg(lx=lgef(x),tx);y[1]=x[1];
  1868.       for(i=2;i<lx;i++) y[i]=(long)centerlift(x[i]);break;
  1869.     case 8: y=cgetg(4,tx);y[1]=copyifstack(x[1]);
  1870.       for(i=2;i<4;i++) y[i]=(long)centerlift(x[i]);break;
  1871.     default: err(lifter1);
  1872.     }
  1873.   return y;
  1874. }
  1875.  
  1876. GEN glt(x,y) GEN x,y; { return gcmp(x,y)<0 ? gun : gzero;}
  1877. GEN gle(x,y) GEN x,y; { return gcmp(x,y)<=0 ? gun : gzero;}
  1878. GEN gge(x,y) GEN x,y; { return gcmp(x,y)>=0 ? gun : gzero;}
  1879. GEN ggt(x,y) GEN x,y; { return gcmp(x,y)>0 ? gun : gzero;}
  1880. GEN geq(x,y) GEN x,y; { return gegal(x,y) ? gun : gzero;}
  1881. GEN gne(x,y) GEN x,y; { return gegal(x,y) ? gzero : gun;}
  1882. GEN gand(x,y) GEN x,y; { return gcmp0(x)||gcmp0(y) ? gzero : gun;}
  1883. GEN gor(x,y) GEN x,y; { return gcmp0(x)&&gcmp0(y) ? gzero : gun;}
  1884.  
  1885. GEN glength(x)
  1886.   GEN x;
  1887. {
  1888.   switch(typ(x))
  1889.   {
  1890.     case 1:
  1891.     case 10: return stoi(lgef(x)-2);
  1892.     case 2: return stoi(lg(x)-2);
  1893.     case 11: return signe(x) ? stoi(lg(x)-2) : gzero;
  1894.     default: return stoi(lg(x)-lontyp[typ(x)]);
  1895.   }
  1896. }
  1897.  
  1898. GEN matsize(x)
  1899.      GEN x;
  1900. {
  1901.   GEN y;
  1902.  
  1903.   switch(typ(x))
  1904.     {
  1905.     case 17: y=cgetg(3,17);y[1]=un;y[2]=lstoi(lg(x)-1);break;
  1906.     case 18: y=cgetg(3,17);y[1]=lstoi(lg(x)-1);y[2]=un;break;
  1907.     case 19: y=cgetg(3,17);y[2]=lstoi(lg(x)-1);
  1908.       y[1]=(lg(x)>1)?lstoi(lg(x[1])-1):zero;break;
  1909.     default: err(matler1);
  1910.     }
  1911.   return y;
  1912. }
  1913.  
  1914. GEN geval(x)
  1915.   GEN x;
  1916. {
  1917.   long av, tetpil, tx = typ(x), lx, i;
  1918.   GEN y, z;
  1919.   if (tx < 9) return gcopy(x);
  1920.   if (tx > 13)
  1921.   {
  1922.     lx = lg(x);
  1923.     y = cgetg(lx, tx);
  1924.     for(i = 1; i < lx; i++) y[i] = (long)geval(x[i]);
  1925.     return y;
  1926.   }
  1927.   switch(tx)
  1928.   {
  1929.     case 9:
  1930.       y=cgetg(3,9);y[1]=(long)geval(x[1]);av=avma;
  1931.       z=geval(x[2]);tetpil=avma;y[2]=lpile(av,tetpil,gmod(z,y[1]));
  1932.       return y;
  1933.     case 10:
  1934.       lx = lgef(x); if (lx == 2) return gzero;
  1935.       y = gzero; z = (GEN)varentries[varn(x)]->value;
  1936.       av = avma;
  1937.       for(i = lx-1; i > 1; i--)
  1938.       {
  1939.         tetpil = avma;
  1940.         y = gadd(geval(x[i]), gmul(z, y));
  1941.       }
  1942.       return gerepile(av, tetpil, y);
  1943.     case 11: err(impl, "evaluation of a power series");
  1944.     case 13: return gdiv(geval(x[1]),geval(x[2]));
  1945.   }
  1946. }
  1947.  
  1948. int isexactzero(g)
  1949.      GEN g;
  1950. {
  1951.   long i;
  1952.   switch (typ(g))
  1953.   {
  1954.     case 1: return !signe(g);
  1955.     case 2:
  1956.     case 7:
  1957.     case 11: return 0;
  1958.     case 3:
  1959.     case 9: return isexactzero(g[2]);
  1960.     case 4:
  1961.     case 5:
  1962.     case 13:
  1963.     case 14: return isexactzero(g[1]);
  1964.     case 6: return isexactzero(g[1])&&isexactzero(g[2]);
  1965.     case 8: return isexactzero(g[2])&&isexactzero(g[3]);
  1966.     case 10: for (i=lgef(g)-1;i>1;i--) if (!isexactzero(g[i])) return 0;
  1967.       return 1;
  1968.     case 17:
  1969.     case 18:
  1970.     case 19: for(i=1;i<lg(g);i++) if(!isexactzero(g[i])) return 0;
  1971.       return 1;
  1972.     default: return 0;
  1973.   }
  1974. }
  1975.  
  1976. GEN simplify(x)
  1977.      GEN x;
  1978. {
  1979.   long tx=typ(x),av,tetpil,i,lx;
  1980.   GEN p1,y;
  1981.  
  1982.   switch(tx)
  1983.     {
  1984.     case 1:
  1985.     case 2:
  1986.     case 3:
  1987.     case 4:
  1988.     case 7:
  1989.     case 15:
  1990.     case 16: return gcopy(x);
  1991.     case 5: return gred(x);
  1992.     case 6:
  1993.     case 8: if(isexactzero(x[2])) return simplify(x[1]);
  1994.       else 
  1995.     {
  1996.       y=cgetg(lg(x),tx);y[1]=(long)simplify(x[1]);
  1997.       y[2]=(long)simplify(x[2]);return y;
  1998.     }
  1999.     case 9: y=cgetg(3,9);y[1]=(long)simplify(x[1]);
  2000.       av=avma;p1=gmod(x[2],y[1]);tetpil=avma;
  2001.       y[2]=lpile(av,tetpil,simplify(p1));return y;
  2002.     case 10: lx=lgef(x);if(lx==2) return gzero;
  2003.       if(lx==3) return simplify(x[2]);
  2004.       y=cgetg(lx,tx);y[1]=x[1];for(i=2;i<lx;i++) y[i]=(long)simplify(x[i]);
  2005.       return y;
  2006.     case 11: if (!signe(x)) return gcopy(x);
  2007.       lx=lg(x);
  2008.       y=cgetg(lx,tx);y[1]=x[1];for(i=2;i<lx;i++) y[i]=(long)simplify(x[i]);
  2009.       return y;
  2010.     case 13: y=cgetg(3,13);y[1]=(long)simplify(x[1]);
  2011.       y[2]=(long)simplify(x[2]);return y;
  2012.     case 14: av=avma;y=gred(x);tetpil=avma;
  2013.       return gerepile(av,tetpil,simplify(y));
  2014.     case 17:
  2015.     case 18:
  2016.     case 19: lx=lg(x);y=cgetg(lx,tx);
  2017.       for(i=1;i<lx;i++) y[i]=(long)simplify(x[i]);
  2018.       return y;
  2019.     default: err(simplifyer1);
  2020.     }
  2021. }
  2022.