home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Professional / OS2PRO194.ISO / os2 / progs / pari / pari_137 / src / gen2.c < prev    next >
Text File  |  1992-05-20  |  52KB  |  2,044 lines
  1. /********************************************************************/
  2. /********************************************************************/
  3. /**                                                                **/
  4. /**                +++++++++++++++++++++++++++++++                 **/
  5. /**                +                             +                 **/
  6. /**                +    OPERATIONS GENERIQUES    +                 **/
  7. /**                +     (deuxieme 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. /*                    OPERATIONS PAR VALEUR                        */
  52. /*                                                                 */
  53. /*      parametres : f pointe sur la fonction appelee;             */
  54. /*                 x,y,... ,z pointent sur des GEN;                */
  55. /*                 le dernier parametre recoit le resultat de      */
  56. /*                 l'operation .                                   */
  57. /*                                                                 */
  58. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  59. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  60.  
  61. void gop0z(f,x)
  62.      
  63.      GEN     (*f)(),x;
  64.      
  65. {
  66.   long    avmacourant;
  67.   GEN   p1;
  68.   
  69.   avmacourant=avma;
  70.   p1=(*f)();gaffect(p1,x);
  71.   avma=avmacourant;
  72. }
  73.  
  74.  
  75. /* operation a un parametre   */
  76.  
  77. void gop1z(f,x,y)
  78.      
  79.      GEN     (*f)(),x,y;
  80.      
  81. {
  82.   long    avmacourant;
  83.   GEN     p1;
  84.   
  85.   avmacourant=avma;p1=(*f)(x);gaffect(p1,y);
  86.   avma=avmacourant;
  87. }
  88.  
  89.  
  90. /* operation a deux parametres */
  91.  
  92. void gop2z(f,x,y,z)
  93.      
  94.      GEN     (*f)(),x,y,z;
  95.      
  96. {
  97.   long    avmacourant;
  98.   GEN     p1;
  99.   
  100.   avmacourant=avma;p1=(*f)(x,y);gaffect(p1,z);
  101.   avma=avmacourant;
  102. }
  103.  
  104. void gops2gsz(f,x,s,z)
  105.      
  106.      GEN     (*f)(),x,z;
  107.      long    s;
  108.      
  109. {
  110.   long    avmacourant;
  111.   GEN     p1;
  112.   
  113.   avmacourant=avma;p1=(*f)(x,s);gaffect(p1,z);
  114.   avma=avmacourant;
  115. }
  116.  
  117.  
  118. void gops2sgz(f,s,y,z)
  119.      
  120.      GEN     (*f)(),y,z;
  121.      long    s;
  122.      
  123. {
  124.   long    avmacourant;
  125.   GEN     p1;
  126.   
  127.   avmacourant=avma;p1=(*f)(s,y);gaffect(p1,z);
  128.   avma=avmacourant;
  129. }
  130.  
  131.  
  132. void gops2ssz(f,s,y,z)
  133.      
  134.      GEN     (*f)(),z;
  135.      long    s,y;
  136.      
  137. {
  138.   long    avmacourant;
  139.   GEN     p1;
  140.   
  141.   avmacourant=avma;p1=(*f)(s,y);gaffect(p1,z);
  142.   avma=avmacourant;
  143. }
  144.  
  145. /* operation a trois parametres */
  146.  
  147. void gop3z(f,x,y,z,t)
  148.      
  149.      GEN     (*f)(),x,y,z,t;
  150.      
  151. {
  152.   long    avmacourant;
  153.   GEN     p1;
  154.   
  155.   avmacourant=avma;p1=(*f)(x,y,z);gaffect(p1,t);
  156.   avma=avmacourant;
  157. }
  158.  
  159.  
  160. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  161. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  162. /*                                                                 */
  163. /*            OPERATIONS AVEC DES ENTIERS COURTS                   */
  164. /*                                                                 */
  165. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  166. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  167.  
  168.  
  169. GEN     gopsg2(f,s,y)
  170.      
  171.      long    s;
  172.      GEN     (*f)(),y;
  173.      
  174. {
  175.   long court[3];
  176.   court[0] = 0x1010003;
  177.   affsi(s,court);return (*f)(court,y);
  178. }
  179.  
  180. GEN     gopgs2(f,y,s)
  181.      
  182.      long    s;
  183.      GEN     (*f)(),y;
  184.      
  185. {
  186.   long court[3];
  187.   court[0] = 0x1010003;
  188.   affsi(s,court);return (*f)(y,court);
  189. }
  190.  
  191. long opgs2(f,y,s)
  192.      
  193.      long    s,(*f)();
  194.      GEN     y;
  195.      
  196. {
  197.   long court[3];
  198.   court[0] = 0x1010003;
  199.   affsi(s,court);return (*f)(y,court);
  200. }
  201.  
  202. void gops1z(f,s,y)
  203.      
  204.      long    s;
  205.      GEN     (*f)(),y;
  206.      
  207. {
  208.   long av=avma; gaffect((*f)(s),y); avma=av;
  209. }
  210.  
  211. void     gopsg2z(f,s,y,z)
  212.      
  213.      long    s;
  214.      GEN     (*f)(),y,z;
  215.      
  216. {
  217.   long av, court[3];
  218.   court[0] = 0x1010003;
  219.   affsi(s,court);av=avma;
  220.   gaffect((*f)(court,y),z);avma=av;
  221. }
  222.  
  223. void     gopgs2z(f,y,s,z)
  224.      
  225.      long    s;
  226.      GEN     (*f)(),y,z;
  227.      
  228. {
  229.   long av, court[3];
  230.   court[0] = 0x1010003;
  231.   affsi(s,court);av=avma;
  232.   gaffect((*f)(y,court),z);avma=av;
  233. }
  234.  
  235. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  236. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  237. /*                                                                 */
  238. /*            CREATION D'UN P-ADIQUE                               */
  239. /*                                                                 */
  240. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  241. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  242.  
  243.  
  244. GEN cgetp(x)
  245.      GEN x;
  246.      
  247. {
  248.   GEN y;
  249.   
  250.   y=cgetg(5,7);
  251.   y[2] = x[2]; y[3] = lcopy(x[3]);
  252.   setprecp(y, precp(x)); y[4] = lgeti(lg(x[3]));
  253.   return y;
  254. }
  255.  
  256. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  257. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  258. /*                                                                 */
  259. /*                       CLONE ET COPIE                            */
  260. /*                                                                 */
  261. /*            Cree la replique exacte d'un GEN existant            */
  262. /*                                                                 */
  263. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  264. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  265.  
  266. GEN     gcopy(x)
  267.      
  268.      GEN   x;
  269. {
  270.   long tx=typ(x),lx=lg(x),i;
  271.   GEN  y;
  272.   
  273.   y=cgetg(lx,tx);
  274.   if(tx<=2)
  275.     for(i=1;i<lx;i++) y[i]=x[i];
  276.   else
  277.   {
  278.     if(tx==10) lx=lgef(x);
  279.     if((tx==11)&&gcmp0(x)) lx=2;
  280.     for(i=1;i<lontyp[tx];i++)
  281.       y[i]=x[i];
  282.     for(i=lontyp[tx];i<lontyp2[tx];i++)
  283.       y[i]=copyifstack(x[i]);
  284.     for(i=lontyp2[tx];i<lx;i++)
  285.       y[i]=lcopy(x[i]);
  286.   }
  287.   return y;
  288. }
  289.  
  290. GEN     forcecopy(x)
  291.      
  292.      GEN   x;
  293. {
  294.   long tx=typ(x),lx=lg(x),i;
  295.   GEN  y;
  296.   
  297.   y=cgetg(lx,tx);
  298.   if(tx<=2)
  299.     for(i=1;i<lx;i++) y[i]=x[i];
  300.   else
  301.   {
  302.     if(tx==10) lx=lgef(x);
  303.     if((tx==11)&&gcmp0(x)) lx=2;
  304.     for(i=1;i<lontyp[tx];i++)
  305.       y[i]=x[i];
  306.     for(i=lontyp[tx];i<lontyp2[tx];i++)
  307.       y[i]=(long)forcecopy(x[i]);
  308.     for(i=lontyp2[tx];i<lx;i++)
  309.       y[i]=(long)forcecopy(x[i]);
  310.   }
  311.   return y;
  312. }
  313.  
  314. long taille(x)
  315.   GEN x;
  316. {
  317.   long i, n, lx = lg(x), tx = typ(x);
  318.   if (tx <= 2) return (tx == 1) ? lgef(x) : lx;
  319.   if ((tx == 11) && gcmp0(x)) return 2;
  320.   if (tx == 10) lx = lgef(x);
  321.   n = lx;
  322.   for(i = lontyp[tx]; i < lx; i++) n += taille(x[i]);
  323.   return n;
  324. }
  325.  
  326. GEN brutcopy(x, y)
  327.   GEN x, y;
  328. {
  329.   long i, lx, tx = typ(x);
  330.   GEN z;
  331.   lx = ((tx == 1) || (tx == 10)) ? lgef(x) : lg(x);
  332.   if (tx <= 2)
  333.     for(i = 0; i < lx; i++) y[i] = x[i];
  334.   else
  335.   {
  336.     if ((tx == 11) && gcmp0(x)) lx = 2;
  337.     z = y + lx;
  338.     for(i = 0; i < lontyp[tx]; i++) y[i] = x[i];
  339.     for(i = lontyp[tx]; i < lx; i++)
  340.     {
  341.       y[i] = (long)brutcopy(x[i], z);
  342.       z += taille(x[i]);
  343.     }
  344.   }
  345.   setlg(y,lx);
  346.   return y;
  347. }
  348.  
  349. GEN gclone(x)
  350.   GEN x;
  351. {
  352.   GEN y = newbloc(taille(x));
  353.   return brutcopy(x, y);
  354. }
  355.  
  356. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  357. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  358. /*                                                                 */
  359. /*                            GREFFE                               */
  360. /*                                                                 */
  361. /*            Greffe d'une serie sur un polynome                   */
  362. /*                                                                 */
  363. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  364. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  365.  
  366.  
  367. GEN     greffe(x,l)
  368.      
  369.      long    l;
  370.      GEN     x;
  371.      
  372. {
  373.   long    i,e,k;
  374.   GEN     y;
  375.   
  376.   if (typ(x)!=10) err(grefer1);
  377.   else
  378.   {
  379.     y=cgetg(l,11);
  380.     if (gcmp0(x))
  381.     {
  382.       y[1]=0x7ffe +l;
  383.       for (i=2;i<l;y[i]=x[2],i++);
  384.     }
  385.     else
  386.     {
  387.       e=gval(x,varn(x));y[1]=0x1008000+e;
  388.       k=lgef(x)-e-1;
  389.       if (k<l)
  390.       {
  391.         for (i=2;i<=k;y[i]=x[i+e],i++);
  392.         for (i=k+1;i<l;y[i]=zero,i++);
  393.       }
  394.       else
  395.         for (i=2;i<l;y[i]=x[i+e],i++);
  396.     }
  397.   }
  398.   setvarn(y,varn(x));return y;
  399. }
  400.  
  401.  
  402. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  403. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  404. /*                                                                 */
  405. /*            CONVERSION GEN-->LONG DU C                           */
  406. /*                                                                 */
  407. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  408. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  409.  
  410. long    gtolong(x)
  411.      
  412.      GEN     x;
  413.      
  414. {
  415.   long    y,tx,av;
  416.   
  417.   tx=typ(x);av=avma;
  418.   
  419.   switch(tx)
  420.     {
  421.     case 1 : y=itos(x);break;
  422.     case 2 : 
  423.     case 4 :
  424.     case 5 : y=itos(ground(x));break;
  425.     case 6 : if (gcmp0(x[2])) y=gtolong(x[1]);
  426.     else err(gtolger);
  427.       break;
  428.     case 8 : if (gcmp0(x[3])) y=gtolong(x[2]);
  429.     else err(gtolger);
  430.       break;
  431.     default: err(gtolger);
  432.   }
  433.   avma=av;
  434.   return y;
  435. }
  436.  
  437.  
  438. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  439. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  440. /*                                                                 */
  441. /*                    COMPARAISON  A ZERO                          */
  442. /*                                                                 */
  443. /*        x pointe sur un GEN;la fonction renvoie 1 si le          */
  444. /*                   GEN est nul,0 sinon .                         */
  445. /*                                                                 */
  446. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  447. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  448.  
  449.  
  450. int gcmp0(x)
  451.      
  452.      GEN     x;
  453.      
  454. {
  455.   long  l,i;
  456.   
  457.   switch(typ(x))
  458.   {
  459.     case 1 :
  460.     case 2 :
  461.     case 10:
  462.     case 11: return !signe(x);
  463.       
  464.     case 3 :
  465.     case 9 : return gcmp0(x[2]);
  466.       
  467.     case 4 :
  468.     case 5 : return !signe(x[1]);
  469.     case 13:
  470.     case 14: return gcmp0(x[1]);
  471.       
  472.     case 6 : return gcmp0(x[1])&&gcmp0(x[2]);
  473.       
  474.     case 8 : return gcmp0(x[2])&&gcmp0(x[3]);
  475.       
  476.     case 7 : return !signe(x[4]);
  477.       
  478.     case 15:
  479.     case 16: 
  480.     case 17:
  481.     case 18:
  482.     case 19: l=lg(x);i=1;
  483.       while ((i<l)&&gcmp0(x[i])) i++;
  484.       return i==l;
  485.       
  486.     default: err(gcmper);
  487.   }
  488. }
  489.  
  490.  
  491. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  492. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  493. /*                                                                 */
  494. /*            COMPARAISONS  A     1  et -1                         */
  495. /*                                                                 */
  496. /*          x pointe sur un GEN;la fonction renvoie 1 si le        */
  497. /*          GEN est l'entier 1 (resp.-1),0 sinon .                 */
  498. /*                                                                 */
  499. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  500. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  501.  
  502.  
  503. int gcmp1(x)
  504.      
  505.      GEN     x;
  506.      
  507. {
  508.   
  509.   long    typy=typ(x),l,y;
  510.   GEN     p1;
  511.   
  512.   switch(typy)
  513.   {
  514.     case 1 : return (lgef(x)==3)&&(signe(x)==1)&&(mant(x,1)==1);
  515.     case 2 : l=avma;p1=subrs(x,1);
  516.       y=(!signe(p1));avma=l;
  517.       return y;
  518.     case 3 :
  519.     case 9 : return gcmp1(x[2]);
  520.     case 4 : return gcmp1(x[1])&&gcmp1(x[2]);
  521.     case 5 : return !(cmpii(x[1],x[2]));
  522.     case 6 : return gcmp1(x[1])&&gcmp0(x[2]);
  523.     case 8 : return gcmp1(x[2])&&gcmp0(x[3]);
  524.     case 7 : return !valp(x)&&gcmp1(x[4]);
  525.     case 10 : return (lgef(x)==3)&&gcmp1(x[2]);
  526.     default: return 0;
  527.   }
  528. }
  529.  
  530.  
  531. int gcmp_1(x)
  532.      
  533.      GEN     x;
  534.      
  535. {
  536.   
  537.   long    typy=typ(x),l,y;
  538.   GEN     p1;
  539.   
  540.   switch(typy)
  541.   {
  542.     case 1 : return (lgef(x)==3)&&(signe(x)== -1)&&(mant(x,1)==1);
  543.     case 2 : l=avma;p1=addsr(1,x);
  544.       y=(!signe(p1));avma=l;
  545.       return y;
  546.     case 3 : l=avma;p1=addsi(1,x[2]);
  547.       y=cmpii(p1,x[1]);avma=l;
  548.       return !y;
  549.     case 4 :
  550.     case 5 : l=avma;p1=addii(x[1],x[2]);
  551.       y=signe(p1);avma=l;
  552.       return !y;
  553.     case 6 : return gcmp_1(x[1])&&gcmp0(x[2]);
  554.     case 8 : return gcmp_1(x[2])&&gcmp0(x[3]);
  555.     case 7 : l=avma;p1=addsi(1,x[4]);
  556.       y=cmpii(p1,x[3]);avma=l;
  557.       return !y;
  558.     case 9 : l=avma;p1=gaddsg(1,x[2]);
  559.       y=(!gegal(p1,x[1]))&&(signe(p1));avma=l;
  560.       return !y;
  561.     case 10 : return (lgef(x)==3)&&gcmp_1(x[2]);
  562.     default: return 0;
  563.     }
  564. }
  565.  
  566. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  567. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  568. /*                                                                 */
  569. /*                       COMPARAISON SIGNEE                        */
  570. /*                                                                 */
  571. /*            rend le signe de x-y si ceci a un sens,              */
  572. /*            sinon erreur.                                        */
  573. /*                                                                 */
  574. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  575. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  576.  
  577.  
  578. int   gcmp(x,y)
  579.      
  580.      GEN     x,y;
  581.      
  582. {
  583.   long  tx,ty,f,avmacourant;
  584.   
  585.   tx=typ(x);ty=typ(y);
  586.   if ((tx>5)||(tx==3)||(ty>5)||(ty==3)) err(gcmper);
  587.   if((tx<=2)&&(ty<=2)) return mpcmp(x,y);
  588.   else
  589.   {
  590.     avmacourant=avma;f=gsigne(gsub(x,y));
  591.     avma=avmacourant;return f;
  592.   }
  593. }
  594.  
  595. int lexcmp(x,y)
  596.      GEN x,y;
  597. {
  598.   long tx=typ(x),ty=typ(y),lx,ly,fl,s,i;
  599.   GEN p1;
  600.  
  601.   if(tx>ty) {p1=x;x=y;y=p1;fl=tx;tx=ty;ty=fl;s= -1;} else s=1;
  602.   ly=lg(y);
  603.   if(tx<17)
  604.     {
  605.       if(ty<17) return s*gcmp(x,y);
  606.       if(ly==1) return s;
  607.       fl=lexcmp(x,y[1]);if(fl) return s*fl;
  608.       return (ly>2)?-s:0;
  609.     }
  610.   lx=lg(x);
  611.   if(ly==1) return (lx==1)?0:s;
  612.   if(lx==1) return -s;
  613.   if((ty==19)&&(tx<19))
  614.     {
  615.       fl=lexcmp(x,y[1]);if(fl) return s*fl;
  616.       return (ly>2)?-s:0;
  617.     }
  618.   if(lx>ly) {p1=x;x=y;y=p1;fl=lx;lx=ly;ly=fl;s= -s;}
  619.   fl=0;for(i=1;(i<lx)&&(!fl);i++) fl=lexcmp(x[i],y[i]);
  620.   if(fl) return s*fl;
  621.   return (ly>lx)?-s:0;
  622. }
  623.  
  624. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  625. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  626. /*                                                                 */
  627. /*                            EGALITE                              */
  628. /*                                                                 */
  629. /*            Renvoie 1 si x=y, 0 sinon                            */
  630. /*                                                                 */
  631. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  632. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  633.  
  634.  
  635. int gegal(x,y)
  636.      
  637.      GEN     x,y;
  638.      
  639. {
  640.   long    avmacourant,f,i,tx=typ(x),ty=typ(y),lx,masq=0xff00ffff;
  641.   GEN     p1;
  642.   
  643.   if(tx!=ty)
  644.   {
  645.     avmacourant=avma;p1=gsub(x,y);f=gcmp0(p1);avma=avmacourant;
  646.   }
  647.   else
  648.   {
  649.     switch(tx)
  650.     {
  651.       case 1: if((x[1]&masq)!=(y[1]&masq)) return 0;
  652.         i=2;lx=lgef(x);while((i<lx)&&(x[i]==y[i])) i++;return i==lx;
  653.       case 10: if(x[1]!=y[1]) return 0;
  654.         i=2;lx=lgef(x);while((i<lx)&&(gegal(x[i],y[i]))) i++;return i==lx;
  655.       case 6: return gegal(x[1],y[1])&&gegal(x[2],y[2]);
  656.       case 3:
  657.       case 9: return (x[1]==y[1])?gegal(x[2],y[2]):gegal(x[1],y[1])&&gegal(x[2],y[2]);
  658.       case 8: return gegal(x[1],y[1])&&gegal(x[2],y[2])&&gegal(x[3],y[3]);
  659.       case 4:
  660.       case 5:
  661.       case 13:
  662.       case 14: avmacourant=avma;f=gegal(gmul(x[1],y[2]),gmul(x[2],y[1]));
  663.     avma=avmacourant;return f;
  664.       case 15: return gegal(x[1],y[1])&&gegal(x[2],y[2])&&gegal(x[3],y[3])&&gegal(x[4],y[4]);
  665.       case 16: return gegal(x[1],y[1])&&gegal(x[2],y[2])&&gegal(x[3],y[3]);
  666.       default: avmacourant=avma;p1=gsub(x,y);f=gcmp0(p1);avma=avmacourant;
  667.     }
  668.   }
  669.   return f;
  670. }
  671.  
  672. int polegal(x,y)
  673.      
  674.      GEN    x,y;
  675.      
  676.      /* a usage interne donc pas de verification de types */
  677.      
  678. {
  679.   long i;
  680.   
  681.   if (x[1] != y[1]) return 0;
  682.   for(i = lgef(x) - 1; i > 1; i--) if (!gegal(x[i],y[i])) return 0;
  683.   return 1;
  684. }
  685.  
  686. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  687. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  688. /*                                                                 */
  689. /*                          VALUATION                              */
  690. /*                                                                 */
  691. /*   renvoie le plus grand exposant de p divisant x quand cela a   */
  692. /*   un sens (erreur pour les types reels, integermod et polymod   */
  693. /*   si le module n'est pas divisible par p, et q-adiques pour     */
  694. /*   q!=p. p doit etre de type entier ou polynome.                 */
  695. /*                                                                 */
  696. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  697. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  698.  
  699. long ggval(x,p)
  700.      
  701.      GEN x,p;
  702.      
  703. {
  704.   long  tx=typ(x),tp=typ(p),lx,i,j,vx,v,av=avma,val,limite=(avma+bot)/2,tetpil;
  705.   GEN  p1,p2,reste;
  706.   
  707.   if(isexactzero(x)) err(gvaler2);
  708.   if((tx<9)&&(tp==10)) return 0;
  709.   switch(tx)
  710.     {
  711.     case 1: if(tp!=1) err(gvaler4);
  712.       val=0;
  713.       while (1)
  714.     {
  715.       p1=dvmdii(x,p,&reste);
  716.       if (signe(reste)) break;
  717.       val++;x=p1;
  718.           if (avma<limite) {tetpil=avma;x=gerepile(av,tetpil,gcopy(x));}
  719.     }
  720.       break;
  721.     case 3: p1=cgeti(lgef(x[1]));
  722.       if((tp!=1)||(!mpdivis(x[1],p,p1))) err(gvaler4);
  723.       p2=cgeti(lgef(x[2]));
  724.       if(mpdivis(x[2],p,p2))
  725.     {val=1;while(mpdivis(p1,p,p1)&&mpdivis(p2,p,p2)) val++;}
  726.       else val=0;
  727.       break;
  728.     case 7: if((tp!=1)||(!gegal(p,x[2]))) err(gvaler4);
  729.       val=valp(x);break;
  730.     case 9: if(tp==1) val=ggval(x[2],p);
  731.       else
  732.     {
  733.       if((tp!=10)||(!poldivis(x[1],p,&p1))) err(gvaler4);
  734.       if(poldivis(x[2],p,&p2))
  735.         {val=1;while(poldivis(p1,p,&p1)&&poldivis(p2,p,&p2)) val++;}
  736.       else val=0;
  737.     }
  738.       break;
  739.     case 10: i=2;while (isexactzero(x[i])) i++;
  740.       if(tp==10)
  741.     {
  742.       v=varn(p);vx=varn(x);if(vx>v) return 0;
  743.       if(vx==v)
  744.         {
  745.           if(ismonome(p)) val=i-2;
  746.           else
  747.         {
  748.           val=0;
  749.           while (1)
  750.             {
  751.               p1=poldivres(x,p,&reste);
  752.               if (signe(reste)) break;
  753.               val++;x=p1;
  754.                       if (avma<limite) {tetpil=avma;x=gerepile(av,tetpil,gcopy(x));}
  755.             }
  756.         }
  757.         }
  758.       else
  759.         {
  760.           val=ggval(x[i],p);
  761.           for(j=i+1;j<lgef(x);j++)
  762.         if(!gcmp0(x[j])) val=min(val,ggval(x[j],p));
  763.         }
  764.     }
  765.       else 
  766.     {
  767.       if(tp!=1) err(gvaler4);
  768.       val=ggval(x[i],p);
  769.       for(j=i+1;j<lgef(x);j++)
  770.         if(!gcmp0(x[j])) val=min(val,ggval(x[j],p));
  771.     }
  772.       break;
  773.       
  774.     case 11: if((tp!=10)&&(tp!=11)&&(tp!=1)) err(gvaler4);
  775.       v=gvar(p);vx=varn(x);if(vx>v) return 0;
  776.       if(vx==v) return (long)(valp(x)/ggval(p,polx[v]));
  777.       val=ggval(x[2],p);
  778.       for(j=3;j<lg(x);j++)
  779.         if(!gcmp0(x[j])) val=min(val,ggval(x[j],p));
  780.       break;
  781.     case 4:
  782.     case 5:
  783.     case 13:
  784.     case 14: val=ggval(x[1],p)-ggval(x[2],p);break;
  785.       
  786.     case 2:
  787.     case 15:
  788.     case 16: err(gvaler4);
  789.     case 6:
  790.     case 8:
  791.     case 17:
  792.     case 18:
  793.     case 19: val=VERYBIGINT;lx=lg(x);
  794.       for(j=1;j<lx;j++) val=min(val,ggval(x[j],p));
  795.       break;
  796.     default: err(gvaler4);
  797.   }
  798.   avma=av;return val;
  799. }
  800.  
  801. long pvaluation(x,p,py)
  802.      
  803.      GEN     x,p,*py;
  804.      
  805. {
  806.   long l,tetpil,val,limite = (bot + avma) / 2;
  807.   GEN  p1,reste;
  808.   
  809.   val=0;*py=cgeti(lgef(x));l=avma;
  810.   while (1)
  811.   {
  812.     p1=dvmdii(x,p,&reste);
  813.     if (signe(reste)) break;
  814.     val++;x=p1;
  815.     if (avma < limite)
  816.       {
  817.         tetpil = avma;
  818.         x = gerepile(l, tetpil, gcopy(x));
  819.       }
  820.   }
  821.   affii(x,*py);avma=l;
  822.   return val;
  823. }
  824.  
  825.  
  826. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  827. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  828. /*                                                                 */
  829. /*                            NEGATION                             */
  830. /*                                                                 */
  831. /*            Cree-x,ou x pointe sur un GEN .                      */
  832. /*                                                                 */
  833. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  834. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  835.  
  836.  
  837. GEN     gneg(x)
  838.      
  839.      GEN     x;
  840.      
  841. {
  842.   long  tx,lx,i,sx;
  843.   GEN   y;
  844.   
  845.   if (gcmp0(x))
  846.   {
  847.     y=gcopy(x);
  848.   }
  849.   else
  850.   {
  851.     tx=typ(x);lx=lg(x);
  852.     
  853.     switch(tx)
  854.     {
  855.       case 1 :
  856.       case 2 : y=gcopy(x);
  857.         sx=signe(x);
  858.         sx= -sx;
  859.         setsigne(y,sx);
  860.         break;
  861.         
  862.       case 3 : y=cgetg(3,tx);y[1]=copyifstack(x[1]);
  863.         if(gcmp0(x[2])) y[2]=zero;
  864.         else
  865.           y[2]=lsubii(y[1],x[2]);
  866.         break;
  867.         
  868.       case 9: y=cgetg(3,tx);y[1]=copyifstack(x[1]);
  869.         y[2]=lneg(x[2]);
  870.         break;
  871.         
  872.       case 4 :
  873.       case 5 :
  874.       case 13:
  875.       case 14: y=cgetg(3,tx);y[1]=lneg(x[1]);y[2]=lcopy(x[2]);
  876.         break;
  877.         
  878.       case 7 : y=cgetp(x);
  879.         if(gcmp0(x[4])) {affsi(0,y[4]);y[1]=x[1];}
  880.         else
  881.           {
  882.             subiiz(y[3],x[4],y[4]);
  883.             setvalp(y,valp(x));
  884.           }
  885.         break;
  886.  
  887.       case 8 : y=cgetg(lx,tx);
  888.         for (i=2;i<lx;i++) y[i]=lneg(x[i]);
  889.         y[1]=copyifstack(x[1]);
  890.         break;
  891.         
  892.       case 6 :
  893.       case 17:
  894.       case 18:
  895.       case 19: y=cgetg(lx,tx);
  896.         for (i=1;i<lx;i++) y[i]=lneg(x[i]);
  897.         break;
  898.  
  899.       case 15:
  900.       case 16: err(gneger);
  901.         
  902.       case 10: lx=lgef(x);
  903.         y=cgetg(lx,tx);
  904.         for (i=2;i<lx;i++) y[i]=lneg(x[i]);
  905.         y[1]=x[1];
  906.         break;
  907.         
  908.       case 11: y=cgetg(lx,tx);
  909.         for (i=2;i<lx;i++) y[i]=lneg(x[i]);
  910.         y[1]=x[1];
  911.     }
  912.   }
  913.   return y;
  914. }
  915.  
  916.  
  917. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  918. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  919. /*                                                                 */
  920. /*                        VALEUR ABSOLUE                           */
  921. /*                                                                 */
  922. /*            Cree un GEN pointant sur la valeur absolue de x      */
  923. /*            a condition que x soit un entier,ou un reel,         */
  924. /*            ou une fraction,ou un complexe;sinon erreur .        */
  925. /*                                                                 */
  926. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  927. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  928.  
  929. GEN     gabs(x,prec)
  930.      GEN     x;
  931.      long   prec;
  932.      
  933. {
  934.   long  tx=typ(x),lx=lg(x),i,l,tetpil;
  935.   GEN   y,p1;
  936.   
  937.   switch(tx)
  938.   {
  939.     case 1 :
  940.     case 2 : y=mpabs(x);break;
  941.       
  942.     case 4 :
  943.     case 5 : y=cgetg(lx,tx);
  944.       y[1]=lmpabs(x[1]);
  945.       y[2]=lmpabs(x[2]);
  946.       break;
  947.       
  948.     case 6 : l=avma;p1=gnorm(x);tetpil=avma;
  949.       y=gsqrt(p1,prec);y=gerepile(l,tetpil,y);
  950.       break;
  951.       
  952.     case 8 : y=transc(gabs,x,prec);break;
  953.       
  954.     case 17:
  955.     case 18:
  956.     case 19: y=cgetg(lx,tx);
  957.       for(i=1;i<lx;i++)
  958.         y[i]=labs(x[i]);
  959.       break;
  960.       
  961.     default: err(gabser);
  962.   }
  963.   return y;
  964. }
  965.  
  966. GEN     gmax(x,y)
  967.      GEN   x,y;
  968.      
  969. {
  970.   return (gcmp(x,y)>=0) ? gcopy(x) : gcopy(y);
  971. }
  972.  
  973. GEN     gmin(x,y)
  974.      GEN   x,y;
  975.      
  976. {
  977.   return (gcmp(x,y)<=0) ? gcopy(x) : gcopy(y);
  978. }
  979.  
  980.  
  981. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  982. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  983. /*                                                                 */
  984. /*                  AFFECTATION  S--> GEN                          */
  985. /*                                                                 */
  986. /*            Etant donnes un long et un GEN,affecte le long       */
  987. /*            dans le GEN,permettant une initialisation .          */
  988. /*                                                                 */
  989. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  990. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  991.  
  992.  
  993. void gaffsg(s,x)
  994.      
  995.      long    s;
  996.      GEN     x;
  997.      
  998. {
  999.   long  tx,l,i,v,val;
  1000.   GEN   reste,p1;
  1001.   
  1002.   tx=typ(x);
  1003.   
  1004.   switch(tx)
  1005.   {
  1006.     case 1 : affsi(s,x);break;
  1007.     case 2 : affsr(s,x);break;
  1008.     case 3 : modsiz(s,x[1],x[2]);break;
  1009.     case 4 :
  1010.     case 5 : affsi(s,x[1]);
  1011.       affsi(1,x[2]);
  1012.       break;
  1013.     case 6 : gaffsg(s,x[1]);
  1014.       gaffsg(0,x[2]);
  1015.       break;
  1016.     case 7 : if (!s) {setvalp(x,defaultpadicprecision);x[4]=zero;}
  1017.     else
  1018.     {
  1019.       val=0;
  1020.       l=avma;
  1021.       while (1)
  1022.       {
  1023.         p1=dvmdsi(s,x[2],&reste);
  1024.         if (signe(reste)) break;
  1025.         val++;s=itos(p1);
  1026.       }
  1027.       avma=l;
  1028.       setvalp(x,val);
  1029.       if (s>0) affsi(s,x[4]);
  1030.       else addsiz(s,x[3],x[4]);
  1031.     }
  1032.     break;
  1033.       
  1034.     case 8 : gaffsg(s,x[2]);
  1035.       gaffsg(0,x[3]);
  1036.       break;
  1037.       
  1038.     case 9 : gaffsg(s,x[2]);break;
  1039.       
  1040.     case 10: v=varn(x);
  1041.       if (!s) x[1]=2;
  1042.       else
  1043.       {
  1044.         x[1]=0x1000003;gaffsg(s,x[2]);
  1045.       }
  1046.       setvarn(x,v);break;
  1047.       
  1048.     case 11: v=varn(x);gaffsg(s,x[2]);l=lg(x);
  1049.       if (!s) x[1]=0x7ffe +l;
  1050.       else
  1051.       {
  1052.         x[1]=0x1008000;
  1053.         for (i=3;i<l;gaffsg(0,x[i]),i++);
  1054.       }
  1055.       setvarn(x,v);break;
  1056.       
  1057.     case 13:
  1058.     case 14: gaffsg(s,x[1]);gaffsg(1,x[2]);
  1059.       break;
  1060.       
  1061.     case 17:
  1062.     case 18: if (lg(x)!=2) err(gaffsger1);
  1063.     else gaffsg(s,x[1]);break;
  1064.       
  1065.     case 19: if (lg(x)!=2) err(gaffsger2);
  1066.     else gaffsg(s,x[1]);break;
  1067.       
  1068.     default: err(gaffer1);
  1069.   }
  1070. }
  1071.  
  1072.  
  1073. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1074. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1075. /*                                                                 */
  1076. /*                    AFFECTATION GENERALE                         */
  1077. /*                                                                 */
  1078. /*            Etant donnes deux GEN x et y,affecte le contenu      */
  1079. /*            de x dans y,si cela est possible .                   */
  1080. /*                                                                 */
  1081. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1082. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1083.  
  1084.  
  1085. void gaffect(x,y)
  1086.      
  1087.      GEN     x,y;
  1088.      
  1089. {
  1090.   long  i,j,k,l,v,vy,val,lx,ly,tx,ty,d,avmacourant;
  1091.   long  num,den;
  1092.   GEN   p1;
  1093.   
  1094.   tx=typ(x);ty=typ(y);lx=lg(x);ly=lg(y);
  1095.   
  1096.   
  1097.   if (tx<10)
  1098.   {
  1099.     if (ty>=10)
  1100.     {
  1101.       switch(ty)
  1102.       {
  1103.         case 10: v=varn(y);gaffect(x,y[2]);
  1104.           for(i=3;i<ly;gaffsg(0,y[i]),i++);
  1105.           if (gcmp0(x)) y[1]=2;
  1106.           else y[1]=0x1000003;
  1107.           setvarn(y,v);break;
  1108.           
  1109.         case 11: v=varn(y);gaffect(x,y[2]);
  1110.           if (gcmp0(x)) y[1]=0x7ffe +ly;
  1111.           else
  1112.           {
  1113.             y[1]=0x1008000;
  1114.             for (i=3;i<ly;gaffsg(0,y[i]),i++);
  1115.           }
  1116.           setvarn(y,v);break;
  1117.           
  1118.         case 13:
  1119.         case 14: gaffect(x,y[1]);gaffsg(1,y[2]); break;
  1120.         case 17:
  1121.         case 18:
  1122.         case 19: err(gaffer2);
  1123.           
  1124.         default: err(gaffer1);
  1125.         }
  1126.     }
  1127.     else
  1128.     {
  1129.       switch(tx)
  1130.       {
  1131.         case 1 :
  1132.           switch(ty)
  1133.           {
  1134.             case 1 : affii(x,y);break;
  1135.             case 2 : affir(x,y);break;
  1136.             case 3 : modiiz(x,y[1],y[2]);break;
  1137.             case 4 :
  1138.             case 5 : affii(x,y[1]);affsi(1,y[2]); break;
  1139.             case 6 : gaffect(x,y[1]);gaffsg(0,y[2]); break;
  1140.             case 7 : if (!signe(x)) {setvalp(y,defaultpadicprecision);y[4]=zero;}
  1141.               else
  1142.               {
  1143.                 l=avma;
  1144.                 val=pvaluation(x,y[2],&p1);
  1145.                 setvalp(y,val);modiiz(p1,y[3],y[4]);
  1146.                 avma=l;
  1147.               }
  1148.               break;
  1149.             case 8 : gaffect(x,y[2]);gaffsg(0,y[3]); break;
  1150.         case 9 : gaffect(x,y[2]);break;
  1151.             default: err(gaffer1);
  1152.             }
  1153.             break;
  1154.           
  1155.         case 2 :
  1156.           switch(ty)
  1157.           {
  1158.             case 1 : err(gaffer3);
  1159.             case 2 : affrr(x,y);break;
  1160.             case 3 :
  1161.             case 4 :
  1162.             case 5 : err(gaffer3);
  1163.             case 6 : gaffect(x,y[1]);gaffsg(0,y[2]); break;
  1164.             case 7 :
  1165.             case 8 : err(gaffer3);
  1166.         case 9 : gaffect(x,y[2]);break;
  1167.             default: err(gaffer1);
  1168.           }
  1169.           break;
  1170.           
  1171.         case 3 :
  1172.           switch(ty)
  1173.           {
  1174.             case 1 :
  1175.             case 2 : err(gaffer4);
  1176.             case 3 : if (!divise(x[1],y[1])) err(gaffer5);
  1177.               modiiz(x[2],y[1],y[2]); break;
  1178.             case 4 :
  1179.             case 5 :
  1180.             case 6 :
  1181.             case 8 : err(gaffer4);
  1182.             case 7 : err(gaffer6);
  1183.         case 9 : gaffect(x,y[2]);break;
  1184.             default: err(gaffer1);
  1185.           }
  1186.           break;
  1187.         case 4 :
  1188.           switch(ty)
  1189.           {
  1190.             case 1 : i=mpdivis(x[1],x[2],y);
  1191.               if (!i) err(gaffer7);
  1192.               break;
  1193.             case 2 : diviiz(x[1],x[2],y);break;
  1194.               
  1195.             case 3 : avmacourant=avma;
  1196.               p1=mpinvmod(x[2],y[1]);
  1197.               modiiz(mulii(x[1],p1),y[1],y[2]);
  1198.               avma=avmacourant;
  1199.               break;
  1200.             case 4 :
  1201.             case 5 : for (i=1;i<=2;affii(x[i],y[i]),i++);break;
  1202.             case 6 : gaffect(x,y[1]);gaffsg(0,y[2]); break;
  1203.             case 7 : if(!signe(x[1]))
  1204.               {
  1205.                 setvalp(y,defaultpadicprecision);y[4]=zero;
  1206.               }
  1207.               else
  1208.               {
  1209.                 l=avma;num=x[1];den=x[2];
  1210.                 if(!(val=pvaluation(num,y[2],&num)))
  1211.                   val= -pvaluation(den,y[2],&den);
  1212.                 p1=mulii(num,mpinvmod(den,y[3]));
  1213.                 modiiz(p1,y[3],y[4]);avma=l;
  1214.                 setvalp(y,val);
  1215.               }
  1216.               break;
  1217.             case 8 : gaffect(x,y[2]);gaffsg(0,y[3]); break;
  1218.         case 9 : gaffect(x,y[2]);break;
  1219.             default: err(gaffer1);
  1220.           }
  1221.           break;
  1222.         case 5 :
  1223.           switch(ty)
  1224.           {
  1225.             case 1 : i=mpdivis(x[1],x[2],y);
  1226.               if (!i) err(gaffer7);
  1227.               break;
  1228.             case 2 : diviiz(x[1],x[2],y);break;
  1229.               
  1230.             case 3 : avmacourant=avma;gaffect(gred(x),y);
  1231.               avma=avmacourant;
  1232.               break;
  1233.             case 4 : gredz(x,y);break;
  1234.               
  1235.             case 5 : for (i=1;i<=2;affii(x[i],y[i]),i++);break;
  1236.             case 6 : gaffect(x,y[1]);gaffsg(0,y[2]); break;
  1237.             case 7 : if(!signe(x[1]))
  1238.               {
  1239.                 setvalp(y,defaultpadicprecision);y[4]=zero;
  1240.               }
  1241.               else
  1242.               {
  1243.                 l=avma;num=x[1];den=x[2];
  1244.                 val=pvaluation(num,y[2],&num)-pvaluation(den,y[2],&den);
  1245.                 p1=mulii(num,mpinvmod(den,y[3]));
  1246.                 modiiz(p1,y[3],y[4]);avma=l;
  1247.                 setvalp(y,val);
  1248.               }
  1249.               break;
  1250.             case 8 : gaffect(x,y[2]);gaffsg(0,y[3]); break;
  1251.         case 9 : gaffect(x,y[2]);break;
  1252.             default: err(gaffer1);
  1253.             }
  1254.             break;
  1255.             
  1256.         case 6 :
  1257.           switch(ty)
  1258.           {
  1259.             case 1 :
  1260.             case 2 :
  1261.             case 3 :
  1262.             case 4 :
  1263.             case 5 :
  1264.             case 7 :
  1265.             case 8 : if (!gcmp0(x[2])) err(gaffer8);
  1266.               gaffect(x[1],y);
  1267.               break;
  1268.               
  1269.             case 6 : for (i=1;i<=2;gaffect(x[i],y[i]),i++);break;
  1270.         case 9 : gaffect(x,y[2]);break;
  1271.             default: err(gaffer1);
  1272.           }
  1273.           break;
  1274.           
  1275.         case 7 :
  1276.           switch(ty)
  1277.           {
  1278.             case 1 :
  1279.             case 2 : err(gaffer10);
  1280.             case 3 : if(valp(x)<0) err(gaffer11);
  1281.               l=avma;
  1282.               val=pvaluation(y[1],x[2],&p1);
  1283.               if(!gcmp1(p1)) err(gaffer11);
  1284.               if(val>(signe(x[4])?(precp(x)+valp(x)):valp(x)+1)) err(gaffer11);
  1285.               modiiz(x[4],y[1],y[2]);
  1286.               avma=l;break;
  1287.             case 4 :
  1288.             case 5 :
  1289.             case 6 :
  1290.             case 8 : err(gaffer10);
  1291.               
  1292.             case 7 : if(cmpii(x[2],y[2])) err(gaffer12);
  1293.               modiiz(x[4],y[3],y[4]);
  1294.               setvalp(y,valp (x));
  1295.               break;
  1296.               
  1297.         case 9 : gaffect(x,y[2]);break;
  1298.             default: err(gaffer1);
  1299.           }
  1300.           break;
  1301.           
  1302.         case 8 :
  1303.           switch(ty)
  1304.           {
  1305.             case 1 :
  1306.             case 3 :
  1307.             case 4 :
  1308.             case 5 :
  1309.             case 7 : if(!gcmp0(x[3])) err(gaffer9);
  1310.               gaffect(x[2],y);
  1311.               break;
  1312.             case 2 : l=avma;p1=co8(x,ly);gaffect(p1,y);
  1313.               avma=l;break;
  1314.             case 6 : ly=precision(y);
  1315.               if(ly)
  1316.               {
  1317.                 l=avma;p1=co8(x,ly);gaffect(p1,y);
  1318.                 avma=l;
  1319.               }
  1320.               else
  1321.               {
  1322.                 if(!gcmp0(x[3])) err(gaffer9);
  1323.                 gaffect(x[2],y);
  1324.               }
  1325.               break;
  1326.             case 8 : if(!gegal(x[1],y[1])) err(gaffer21);
  1327.               for (i=2;i<=3;gaffect(x[i],y[i]),i++);
  1328.               break;
  1329.         case 9 : gaffect(x,y[2]);break;
  1330.             default: err(gaffer1);
  1331.           }
  1332.           break;
  1333.     case 9 : if(ty!=9) err(gaffer17);
  1334.       if (gdivise(x[1],y[1]))
  1335.         gmodz(x[2],y[1],y[2]);else err(gaffer18);
  1336.       break;
  1337.         default: err(gaffer1);
  1338.     }
  1339.     }
  1340.   }
  1341.   else
  1342.   {
  1343.     if (ty<9) err(gaffer13);
  1344.     switch(tx)
  1345.     {
  1346.       case 10: v=varn(x);
  1347.         switch(ty)
  1348.         {
  1349.           case 10: if((vy=varn(y))>v) err(gaffer13);
  1350.             if(vy==v)
  1351.             {
  1352.               d=lgef(x);
  1353.               if (d>ly) err(gaffer14);
  1354.               y[1]=x[1];
  1355.               for (i=2;i<d;gaffect(x[i],y[i]),i++);
  1356.             }
  1357.             else
  1358.             {
  1359.               gaffect(x,y[2]);
  1360.               for(i=3;i<ly;gaffsg(0,y[i]),i++);
  1361.               if (signe(x)) y[1]=0x1000003;
  1362.               else y[1]=2;
  1363.               setvarn(y,vy);
  1364.             }
  1365.             break;
  1366.             
  1367.           case 11: if((vy=varn(y))>v) err(gaffer13);
  1368.             if (!signe(x)) gaffsg(0,y);
  1369.             else
  1370.             {
  1371.               if(vy==v)
  1372.               {
  1373.                 i=gval(x,v);setvalp(y,i);setsigne(y,1);
  1374.                 k=lgef(x)-i;
  1375.                 if(k>ly) k=ly;
  1376.                 for (j=2;j<k;j++)
  1377.                   gaffect(x[i+j],y[j]);
  1378.                 for (j=k;j<ly;j++)
  1379.                   gaffsg(0,y[j]);
  1380.               }
  1381.               else
  1382.               {
  1383.                 gaffect(x,y[2]);
  1384.                 if (!signe(x)) y[1]=0x7ffe +ly;
  1385.                 else
  1386.                 {
  1387.                   y[1]=0x1008000;
  1388.                   for (i=3;i<ly;gaffsg(0,y[i]),i++);
  1389.                 }
  1390.                 setvarn(y,vy);
  1391.               }
  1392.             }
  1393.             break;
  1394.             
  1395.           case 9: gmodz(x,y[1],y[2]);break;
  1396.           case 13:
  1397.           case 14: gaffect(x,y[1]);gaffsg(1,y[2]); break;
  1398.           case 17:
  1399.           case 18:
  1400.           case 19: err(gaffer15);
  1401.             
  1402.           default: err(gaffer1);
  1403.         }
  1404.         break;
  1405.         
  1406.       case 11: if (ty!=11) err(gaffer16);
  1407.         v=varn(x);if((vy=varn(y))>v) err(gaffer13);
  1408.         if(vy==v)
  1409.         {
  1410.           y[1]=x[1];
  1411.           if(!gcmp0(x))
  1412.           {
  1413.             k=lx;if (k>ly) k=ly;
  1414.             for (i=2;i<k;i++) gaffect(x[i],y[i]);
  1415.             for (i=k;i<ly;i++) gaffsg(0,y[i]);
  1416.           }
  1417.         }
  1418.         else
  1419.         {
  1420.           gaffect(x,y[2]);
  1421.           if (!signe(x)) y[1]=0x7ffe +ly;
  1422.           else
  1423.           {
  1424.             y[1]=0x1008000;
  1425.             for (i=3;i<ly;gaffsg(0,y[i]),i++);
  1426.           }
  1427.           setvarn(y,vy);
  1428.         }
  1429.         break;
  1430.         
  1431.       case 13:
  1432.       case 14: switch(ty)
  1433.       {
  1434.         case 10:
  1435.         case 17:
  1436.         case 18:
  1437.         case 19: err(gaffer19);
  1438.         case 11: gdivz(x[1],x[2],y);break;
  1439.         case 9: avmacourant=avma;p1=ginvmod(x[2],y[1]);
  1440.           gmod(gmul(x[1],p1),y[1],y[2]);
  1441.           avma=avmacourant;
  1442.           break;
  1443.           
  1444.         case 13:
  1445.           if (tx==14) gredz(x,y);
  1446.           else
  1447.             for (i=1;i<=2;gaffect(x[i],y[i]),i++);
  1448.           break;
  1449.           
  1450.         case 14: for (i=1;i<=2;gaffect(x[i],y[i]),i++); break;
  1451.           
  1452.         default: err(gaffer1);
  1453.       }
  1454.       break;
  1455.  
  1456.       case 15:
  1457.       case 16:
  1458.       case 17:
  1459.       case 18:
  1460.       case 19: if (ty!=tx) err(gaffer20);
  1461.         if (ly!=lx) err(gaffer20);
  1462.         for (i=1;i<lx;gaffect(x[i],y[i]),i++);
  1463.         break;
  1464.         
  1465.       default: err(gaffer1);
  1466.     }
  1467.   }
  1468. }
  1469.  
  1470. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1471. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1472. /*                                                                 */
  1473. /*                 CONVERSION QUAD-->REEL,COMPLEXE                 */
  1474. /*                          OU P-ADIQUE                            */
  1475. /*                                                                 */
  1476. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1477. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1478.  
  1479. GEN     co8(x,l)
  1480.      
  1481.      GEN     x;
  1482.      long    l;
  1483.      
  1484. {
  1485.   long  av,tetpil;
  1486.   GEN   y,p1,p2,p3;
  1487.   
  1488.   p3=(GEN)x[1];
  1489.   av=avma;p1=gmul2n(p3[3],-1);
  1490.   p2=gsqrt(gsub(gmul(p1,p1),p3[2]),l);
  1491.   p1=gmul(x[3],gsub(p2,p1));
  1492.   tetpil=avma;y=gerepile(av,tetpil,gadd(p1,x[2]));
  1493.   return y;
  1494. }
  1495.  
  1496. GEN   cvtop(x,p,l)
  1497.      
  1498.      GEN  x,p;
  1499.      long l;
  1500.      
  1501. {
  1502.   GEN z,p1,p2,p3;
  1503.   long av,tetpil,n,tx=typ(x);
  1504.   
  1505.   if(typ(p)!=1) err(cvtoper1);
  1506.   if(gcmp0(x)) z=ggrando(p,l);
  1507.   else
  1508.   {
  1509.     switch(tx)
  1510.     {
  1511.       case 1: z=gadd(x,ggrando(p,ggval(x,p)+l));break;
  1512.       case 2: err(cvtoper2);
  1513.       case 3: n=ggval(x[1],p);
  1514.         if(n>l) n=l;
  1515.         z=gadd(x[2],ggrando(p,n));break;
  1516.       case 4:
  1517.       case 5: n=ggval(x[1],p);
  1518.         n-=ggval(x[2],p);
  1519.         z=gadd(x,ggrando(p,n+l));break;
  1520.       case 6: av=avma;p1=gsqrt(gaddgs(ggrando(p,l),-1));
  1521.         p1=gmul(p1,x[2]);tetpil=avma;
  1522.         z=gerepile(av,tetpil,gadd(p1,x[1]));
  1523.         break;
  1524.       case 7: z=gprec(x,l);break;
  1525.       case 8:
  1526.         av=avma;p1=(GEN)x[1];p3=gmul2n(p1[3],-1);
  1527.         p2=gsub(gmul(p3,p3),p1[2]);
  1528.         switch(typ(p2))
  1529.           {
  1530.           case 1: n=ggval(p2,p);
  1531.             p2=gadd(p2,ggrando(p,n+l));break;
  1532.           case 3: break;
  1533.           case 4:
  1534.           case 5: n=ggval(p2[1],p);
  1535.             n-=ggval(p2[2],p);
  1536.             p2=gadd(p2,ggrando(p,n+l));break;
  1537.           case 7: break;
  1538.           default: err(gadder14);
  1539.           }
  1540.         p2=gsqrt(p2);
  1541.         p1=gmul(x[3],gsub(p2,p3));tetpil=avma;
  1542.         z=gerepile(av,tetpil,gadd(x[2],p1));
  1543.         break;
  1544.       default: err(cvtoper2);
  1545.     }
  1546.   }
  1547.   return z;
  1548. }
  1549.  
  1550. GEN gcvtop(x,p,r)
  1551.      GEN x,p;
  1552.      long r;
  1553.  
  1554. {
  1555.   long i,tx=typ(x),lx;
  1556.   GEN y;
  1557.  
  1558.   if(tx<9) return cvtop(x,p,r);
  1559.   switch(tx)
  1560.     {
  1561.     case 10: lx=lgef(x);y=cgetg(lx,10);y[1]=x[1];
  1562.       for(i=2;i<lx;i++) y[i]=(long)gcvtop(x[i],p,r);break;
  1563.     case 11:
  1564.       if(gcmp0(x)) y=gcopy(x);
  1565.       else 
  1566.     {
  1567.       lx=lg(x);y=cgetg(lx,11);y[1]=x[1];
  1568.       for(i=2;i<lx;i++) y[i]=(long)gcvtop(x[i],p,r);
  1569.     }
  1570.       break;
  1571.     case 9:
  1572.     case 13:
  1573.     case 14:
  1574.     case 17:
  1575.     case 18:
  1576.     case 19: lx=lg(x);y=cgetg(lx,tx);
  1577.       for(i=1;i<lx;i++) y[i]=(long)gcvtop(x[i]);break;
  1578.     default: err(cvtoper2);
  1579.     }
  1580.   return y;
  1581. }
  1582.  
  1583. long   gexpo(x)
  1584.      GEN   x;
  1585.      
  1586. {
  1587.   long tx=typ(x),lx=lg(x),e,i,y,av;
  1588.   GEN  p1;
  1589.   
  1590.   switch(tx)
  1591.   {
  1592.     case 1 :if(!signe(x)) err(gexpoer2);
  1593.       return expi(x);
  1594.     case 4 :
  1595.     case 5 : if(!signe(x[1])) err(gexpoer2);
  1596.       return expi(x[1])-expi(x[2]);
  1597.     case 2 : return expo(x);
  1598.     case 6 : av=avma;p1=gnorm(x);y=(gexpo(p1)>>1); avma=av;break;
  1599.     case 8 :  if(gcmp0(x)) err(gexpoer2);
  1600.       av=avma;p1=cgetg(3,6);p1[1]=lgetr(3);
  1601.       p1[2]=lgetr(3);gaffect(x,p1);y=gexpo(p1);
  1602.       avma=av;break;
  1603.     case 10: lx=lgef(x);
  1604.     case 11: 
  1605.     case 17:
  1606.     case 18:
  1607.     case 19: y= -BIGINT;
  1608.       for(i=lontyp[tx];i<lx;i++)
  1609.       {
  1610.         e=gexpo(x[i]);if(e>y) y=e;
  1611.       }
  1612.       break;
  1613.     default: err(gexpoer1);
  1614.   }
  1615.   return y;
  1616. }
  1617.  
  1618. long gsize(x)
  1619.      GEN x;
  1620. {
  1621.   long l;
  1622.   if(gcmp0(x)) return 0;
  1623.   l=gexpo(x)*L2SL10;
  1624.   return l;
  1625. }
  1626.  
  1627. void normalize(px)
  1628.      
  1629.      GEN  *px;
  1630.      
  1631. {
  1632.   long    i,j,v,l,tetpil,e,lx,f;
  1633.   GEN     p1,x;
  1634.   
  1635.   if(typ(x= *px)!=11) err(gnormaler);
  1636.   i=2;f=1;lx=lg(x);v=varn(x);
  1637.   while (f&&(i<lx))
  1638.     {f=isexactzero(x[i]);i++;}
  1639. /*    {f=gcmp0(x[i]);i++;} */
  1640.   i--;
  1641.   if(i>2)
  1642.   {
  1643.     l=(long)(x+lx);e=valp(x);
  1644.     if (f)
  1645.     {
  1646.       avma=l;
  1647.       p1=cgetg(3,11);p1[1]=0x7ffe +lx;
  1648.       *px=p1;
  1649.     }
  1650.     else
  1651.     {
  1652.       tetpil=avma;
  1653.       p1=cgetg(lx-i+2,11);
  1654.       p1[1]=0x1007ffe +e+i;
  1655.       for (j=i;j<lx;j++)
  1656.         p1[j-i+2]=lcopy(x[j]);
  1657.       *px=gerepile(l,tetpil,p1);
  1658.     }
  1659.     setvarn(*px,v);
  1660.   }
  1661.   else
  1662.   {
  1663.     if(f) setsigne(x,0);else setsigne(x,1);
  1664.   }
  1665. }
  1666.  
  1667. void normalizepol(px)
  1668.      
  1669.      GEN  *px;
  1670.      
  1671. {
  1672.   long    i,lx,f;
  1673.   GEN     x;
  1674.   
  1675.   if(typ(x= *px)!=10) err(gnormaler);
  1676.   lx=lgef(x);
  1677.   if(lx>2)
  1678.     {
  1679.       f=1;i=lx-1;
  1680.       while (f&&(i>1))
  1681.     {f=isexactzero(x[i]);i--;}
  1682.       if(f) {setlgef(*px,2);setsigne(*px,0);}
  1683.       else
  1684.     {
  1685.       i++;f=gcmp0(x[i]);setlgef(*px,i+1);
  1686.       while(f&&(i>1))
  1687.         {f=gcmp0(x[i]);i--;}
  1688.       setsigne(*px,f?0:1);
  1689.     }
  1690.     }
  1691.   else setsigne(*px,0);
  1692. }
  1693.  
  1694. int gsigne(x)
  1695.      GEN x;
  1696. {
  1697.   switch(typ(x))
  1698.   {
  1699.     case 1:
  1700.     case 2: return signe(x);
  1701.     case 4:
  1702.     case 5: return (signe(x[2])>0) ? signe(x[1]) : -signe(x[1]);
  1703.     default: err(gsigner);
  1704.   }
  1705. }
  1706.  
  1707. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1708. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1709. /*                                                                 */
  1710. /*                            PUISSANCE                            */
  1711. /*                                                                 */
  1712. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1713. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1714.  
  1715. GEN     gpui(x,n,prec)
  1716.      
  1717.      GEN     x,n;
  1718.      long    prec;
  1719.      
  1720. {
  1721.   long    av1,av2,av3,lx,tx,i,j,tetpil;
  1722.   unsigned long m;
  1723.   GEN     p1,p2,y,z,alp;
  1724.   
  1725.   if (typ(n)==1)
  1726.   {
  1727.     if((lgef(n)<=3)&&cmpis(n,0x7fffffff)<0) y=gpuigs(x,itos(n),prec);
  1728.     else
  1729.     {
  1730.       z=x;av1=avma;
  1731.       if((typ(x)!=16)&&(typ(x)!=15)) y=gun;
  1732.       else
  1733.       {
  1734.     if(typ(x)==16)
  1735.       {
  1736.         p1=mulii(x[1],x[3]);p2=shifti(mulii(x[2],x[2]),-2);
  1737.         y=cgetg(4,16);y[1]=un;y[2]=mpodd(x[2]) ? un : zero;
  1738.         y[3]=lsubii(p1,p2);
  1739.       }
  1740.     else
  1741.       {
  1742.         y=cgetg(5,15);y[1]=un;
  1743.         p1=subii(mulii(x[2],x[2]),shifti(mulii(x[1],x[3]),2));
  1744.         p2=racine(p1);if(mpodd(subii(p2,x[2]))) p2=addsi(-1,p2);
  1745.         y[2]=(long)p2;y[3]=lshifti(subii(mulii(p2,p2),p1),-2);
  1746.         p1=cgetr(prec);y[4]=(long)p1;p1[1]=0x800000-((prec-2)<<5);
  1747.         p1[2]=0;tetpil=avma;y=gerepile(av1,tetpil,gcopy(y));
  1748.       }
  1749.       }
  1750.       for (i=lgef(n)-1;i>2;i--)
  1751.       {
  1752.         for (m=n[i],j=0;j<32;j++,m>>=1)
  1753.         {
  1754.           if (m&1) y=gmul(y,z);
  1755.           z=gsqr(z);
  1756.         }
  1757.       }
  1758.       for (m=n[2];m>1;m>>=1)
  1759.       {
  1760.         if (m&1) y=gmul(y,z);
  1761.         z=gsqr(z);
  1762.       }
  1763.       if (signe(n)>0)
  1764.       {
  1765.         tetpil=avma;y=gmul(y,z);
  1766.       }
  1767.       else
  1768.       {
  1769.         y=gmul(y,z);
  1770.         tetpil=avma;y=ginv(y);
  1771.       }
  1772.       y=gerepile(av1,tetpil,y);
  1773.     }
  1774.   }
  1775.   else
  1776.   {
  1777.     if((tx=typ(x))>=17)
  1778.       {
  1779.     y=cgetg(lx=lg(x),tx);for(i=1;i<lx;i++) y[i]=lpui(x[i],n,prec);
  1780.       }
  1781.     else
  1782.       {
  1783.     if(tx!=11)
  1784.       {
  1785.         av1=avma;
  1786.         if(gcmp0(x))
  1787.           {
  1788.         if(gcmpgs(p1=greal(n),0)<=0) err(gpuier3);
  1789.         av2=precision(x);
  1790.         if(av2)
  1791.           {
  1792.             p1=ground(gmulsg(gexpo(x),p1));
  1793.             if(lgef(p1)>3) err(gpuier4);
  1794.             av2=itos(p1);if(abs(av2)>=0x800000) err(gpuier4);
  1795.             avma=av1;y=cgetr(3);y[2]=0;y[1]=0x800000+av2;
  1796.           }
  1797.         else {avma=av1;y=gcopy(x);}
  1798.           }
  1799.         else
  1800.           {
  1801.         y=gmul(n,glog(x,prec));av2=avma;
  1802.         y=gerepile(av1,av2,gexp(y,prec));
  1803.           }
  1804.       }
  1805.     else
  1806.       {
  1807.         if(valp(x)) err(gpuier2);
  1808.         if(gvar9(n)>varn(x))
  1809.           {
  1810.         if(gcmp1(x[2]))
  1811.           {
  1812.             av1=avma;alp=gaddgs(n,1);
  1813.             av2=avma;y=cgetg(lx=lg(x),11);y[1]=0x1008000;
  1814.             y[2]=un;setvarn(y,varn(x));
  1815.             for(i=3;i<lx;i++)
  1816.               {
  1817.             av3=avma;p1=gzero;
  1818.             for(j=1;j<i-1;j++)
  1819.               {
  1820.                 p2=gsubgs(gmulgs(alp,j),i-2);
  1821.                 p1=gadd(p1,gmul(gmul(p2,x[j+2]),y[i-j]));
  1822.               }
  1823.             tetpil=avma;
  1824.             y[i]=lpile(av3,tetpil,gdivgs(p1,i-2));
  1825.               }
  1826.             y=gerepile(av1,av2,y);
  1827.           }
  1828.         else
  1829.           {
  1830.             av1=avma;p1=gdiv(x,x[2]);
  1831.             p1=gpui(p1,n,prec);p2=gpui(x[2],n,prec);
  1832.             tetpil=avma;y=gerepile(av1,tetpil,gmul(p1,p2));
  1833.           }
  1834.           }
  1835.         else
  1836.           {
  1837.         av1=avma;y=gmul(n,glog(x,prec));av2=avma;
  1838.         y=gerepile(av1,av2,gexp(y,prec));
  1839.           }
  1840.       }
  1841.       }
  1842.   }
  1843.   return y;
  1844. }
  1845.  
  1846. GEN     gpuigs(x,n,prec)
  1847.      
  1848.      GEN     x;
  1849.      long    n,prec;
  1850.      
  1851. {
  1852.   long    lx,av,m,dd,tetpil,i,j;
  1853.   GEN     y,z,p1,p2;
  1854.   
  1855.   
  1856.   if (!n)
  1857.   {
  1858.     lx=lg(x);
  1859.     switch(typ(x))
  1860.     {
  1861.       case 1 :
  1862.       case 2 :
  1863.       case 4 :
  1864.       case 5 : y=gun;break;
  1865.       case 6 : y=cgetg(3,6);y[1]=un;
  1866.         y[2]=zero;
  1867.         break;
  1868.       case 3 : y=cgetg(3,3);y[1]=copyifstack(x[1]);
  1869.         y[2]=un;
  1870.         break;
  1871.       case 8 : y=cgetg(4,8);y[1]=copyifstack(x[1]);
  1872.         y[2]=un;y[3]=zero;
  1873.         break;
  1874.       case 10: 
  1875.       case 11: 
  1876.       case 13: 
  1877.       case 14: y=polun[varn(x)]; break;
  1878.       case 9: y=cgetg(3,9);y[1]=copyifstack(x[1]);
  1879.         y[2]=lpuigs(x[2],0);
  1880.         break;
  1881.       case 15 : av=avma;y=cgetg(5,15);y[1]=un;
  1882.         p1=subii(mulii(x[2],x[2]),shifti(mulii(x[1],x[3]),2));
  1883.     p2=racine(p1);if(mpodd(subii(p2,x[2]))) p2=addsi(-1,p2);
  1884.     y[2]=(long)p2;y[3]=lshifti(subii(mulii(p2,p2),p1),-2);
  1885.         p1=cgetr(prec);y[4]=(long)p1;p1[1]=0x800000-((prec-2)<<5);
  1886.     p1[2]=0;tetpil=avma;y=gerepile(av,tetpil,gcopy(y));break;
  1887.       case 16 : y=cgetg(4,16);y[1]=un;y[2]=mpodd(x[2]) ? un : zero;
  1888.         av=avma;p1=mulii(x[1],x[3]);p2=shifti(mulii(x[2],x[2]),-2);tetpil=avma;
  1889.         y[3]=lpile(av,tetpil,subii(p1,p2));break;
  1890.       case 17:
  1891.       case 18: err(gpuier1);
  1892.       case 19: if (lx!=(lg(x[1]))) err(gpuier1);
  1893.         else
  1894.         {
  1895.           y=cgetg(lx,19);
  1896.           for (j=1;j<lx;j++)
  1897.           {
  1898.             y[j]=lgetg(lx,18);
  1899.             for (i=1;i<lx;i++)
  1900.               coeff(y,i,j)=(i!=j) ? zero :
  1901.                 lpuigs(coeff(x,i,i),0);
  1902.           }
  1903.         }
  1904.     }
  1905.   }
  1906.   else if (n==1) y=gcopy(x);
  1907.   else if (n== -1) y=ginv(x);
  1908.   else
  1909.   {
  1910.     m=abs(n);
  1911.     if(ismonome(x))
  1912.     {
  1913.       j=lgef(x);
  1914.       dd=(j-3)*m+3;
  1915.       av=avma;p1=cgetg(dd,10);p1[1]=0x1000000+dd;setvarn(p1,varn(x));
  1916.       p1[dd-1]=lpuigs(x[j-1],m);
  1917.       for(i=2;i<dd-1;i++) p1[i]=zero;
  1918.       if(n<0)
  1919.       {
  1920.         tetpil=avma;y=cgetg(3,13);y[1]=(long)denom(p1[dd-1]);
  1921.         y[2]=lmul(p1,y[1]);y=gerepile(av,tetpil,y);
  1922.       }
  1923.       else y=p1;
  1924.     }
  1925.     else
  1926.     {
  1927.       z=x;av=avma;
  1928.       if((typ(x)!=16)&&(typ(x)!=15)) y=gun;
  1929.       else
  1930.       {
  1931.     if(typ(x)==16)
  1932.       {
  1933.         p1=mulii(x[1],x[3]);p2=shifti(mulii(x[2],x[2]),-2);
  1934.         y=cgetg(4,16);y[1]=un;if(mpodd(x[2])) y[2]=un;else y[2]=zero;
  1935.         y[3]=lsubii(p1,p2);
  1936.       }
  1937.     else
  1938.       {
  1939.         y=cgetg(5,15);y[1]=un;
  1940.         p1=subii(mulii(x[2],x[2]),shifti(mulii(x[1],x[3]),2));
  1941.         p2=racine(p1);if(mpodd(subii(p2,x[2]))) p2=addsi(-1,p2);
  1942.         y[2]=(long)p2;y[3]=lshifti(subii(mulii(p2,p2),p1),-2);
  1943.         p1=cgetr(prec);y[4]=(long)p1;p1[1]=0x800000-((prec-2)<<5);
  1944.         p1[2]=0;tetpil=avma;y=gerepile(av,tetpil,gcopy(y));
  1945.       }
  1946.       }
  1947.       for(;m>1;m>>=1)
  1948.       {
  1949.         if (m&1) y=gmul(y,z);
  1950.         z=gsqr(z);
  1951.       }
  1952.       if (n>0)
  1953.       {
  1954.         tetpil=avma;y=gmul(y,z);
  1955.       }
  1956.       else
  1957.       {
  1958.         y=gmul(y,z);
  1959.         tetpil=avma;y=ginv(y);
  1960.       }
  1961.       y=gerepile(av,tetpil,y);
  1962.     }
  1963.   }
  1964.   return y;
  1965. }
  1966.  
  1967. int iscomplex(x)
  1968.   GEN x;
  1969.   
  1970. {
  1971.   long tx=typ(x);
  1972.   
  1973.   switch(tx)
  1974.   {
  1975.     case 1:
  1976.     case 2:
  1977.     case 4:
  1978.     case 5: return 0;
  1979.     case 6: return !gcmp0(x[2]);
  1980.     case 8: return signe(((GEN)x[1])[2])>0;
  1981.     default: err(iscomplexer1);
  1982.   }
  1983. }
  1984.  
  1985. int ismonome(x)
  1986.      GEN x;
  1987.  
  1988. {
  1989.   long i,l;
  1990.  
  1991.   if((typ(x)!=10)||(!signe(x))) return 0;
  1992.   l=lgef(x)-1;i=2;
  1993.   while((i<=l)&&isexactzero(x[i])) i++;
  1994.   return i==l;
  1995. }
  1996.  
  1997. GEN     gsqr(x)
  1998.      
  1999.      GEN     x;
  2000.      
  2001. {
  2002.   long  tx=typ(x),lx,av,i,j;
  2003.   GEN   y,p1;
  2004.   
  2005.   if(tx==7)
  2006.   {
  2007.     y=cgetg(5,7);
  2008.     y[2]=x[2];
  2009.     if(!cmpis(x[2] ,2))
  2010.     {
  2011.       i=precp(x)+1;av=avma;
  2012.       p1=addii(x[3],shifti(x[4],1));
  2013.       if(!gcmp0(p1))
  2014.       {
  2015.         j=vali(p1);if(j<i) i=j;
  2016.       }
  2017.       avma=av;y[3]=lshifti(x[3],i);
  2018.       setprecp(y,precp(x)+i);
  2019.     }
  2020.     else
  2021.     {
  2022.       y[3]=lcopy(x[3]);
  2023.       setprecp(y,precp(x));
  2024.     }
  2025.     y[4]=lgeti(lg(y[3]));
  2026.     setvalp(y,2*valp(x));av=avma;
  2027.     modiiz(mulii(x[4],x[4]),y[3],y[4]);
  2028.     avma=av;
  2029.   }
  2030.   else if(tx==15) y=sqcompreal(x);
  2031.   else if(tx==16) y=sqcomp(x);
  2032.   else
  2033.   {
  2034.     if((tx<17)||((tx==19)&&(lg(x)==lg(x[1])))) y=gmul(x,x);
  2035.     else
  2036.     {
  2037.       y=cgetg(lx=lg(x),tx);
  2038.       for(i=1;i<lx;i++)
  2039.         y[i]=lsqr(x[i]);
  2040.     }
  2041.   }
  2042.   return y;
  2043. }
  2044.