home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 2 / RISC_DISC_2.iso / pd_share / program / code / pari / c / gen2 < prev    next >
Encoding:
Text File  |  1991-11-28  |  51.8 KB  |  2,121 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. /*                                                                 */
  406. /*            CONVERSION      GEN <-->REEL DU C                    */
  407. /*                                                                 */
  408. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  409. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  410.  
  411.  
  412. double rtodbl(x)
  413.      GEN x;
  414. {
  415.   double y,ma,co=32.0;
  416.   long ex,s=signe(x);
  417.   
  418.   if((!s)||((ex=expo(x))< -1023)) return 0.0;
  419.   if(ex>=0x3ff) err(rtodber);
  420.   ma=(ulong)x[2]+(ulong)x[3]/exp2(co);
  421.   y=exp2(ex+1-co);
  422.   return (s<0) ? -y*ma : y*ma;
  423. }
  424.  
  425.  
  426. GEN dbltor(x)
  427.      double x;
  428. {
  429.   GEN z;
  430.   double co=32.0;
  431.   ulong y;
  432.   long ex;
  433.   int n;
  434.   
  435.   if(x==0) {z=cgetr(3);z[2]=0;z[1]=0x800000-308;return z;}
  436.   z=cgetr(4);if(x<0) {setsigne(z,-1);x= -x;} else setsigne(z,1);
  437.   ex=floor(log2(x));setexpo(z,ex);x=x*exp2(-ex+co-1);y=x;x-=y;z[2]=y;
  438.   y=exp2(co)*x;z[3]=y;
  439.  
  440.   if(!(z[2]&(1<<31)))
  441.     {
  442.       n=bfffo(z[2]);z[2]=shiftl(z[2],n);z[3]=shiftl(z[3],n);
  443.       z[2]+=hiremainder;setexpo(z,expo(z)-n);
  444.     }
  445.   return z;
  446. }
  447.  
  448. double  gtodouble(x)
  449.      
  450.      GEN     x;
  451.      
  452. {
  453.   GEN     x1;
  454.   long  t=typ(x);
  455.   static long reel4[4]={0x2010004,0,0,0};
  456.   
  457.   if (t==2) x1=x;
  458.   else gaffect(x,x1=(GEN)reel4);
  459.   return rtodbl(x1);
  460. }
  461.  
  462. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  463. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  464. /*                                                                 */
  465. /*            CONVERSION GEN-->LONG DU C                           */
  466. /*                                                                 */
  467. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  468. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  469.  
  470. long    gtolong(x)
  471.      
  472.      GEN     x;
  473.      
  474. {
  475.   long    y,tx,av;
  476.   GEN     p1,p2;
  477.   
  478.   tx=typ(x);av=avma;
  479.   
  480.   switch(tx)
  481.     {
  482.     case 1 : y=itos(x);break;
  483.       
  484.     case 4 :
  485.       
  486.     case 5 : p1=dvmdii(x[1],x[2],&p2);
  487.       if (!cmpis(p2,0))
  488.     {
  489.       y=signe(x)*itos(p1);break;
  490.     }
  491.       else err(gtolger);
  492.       
  493.     case 6 : if (gcmp0(x[2])) y=gtolong(x[1]);
  494.     else err(gtolger);
  495.       break;
  496.       
  497.     case 8 : if (gcmp0(x[3])) y=gtolong(x[2]);
  498.     else err(gtolger);
  499.       break;
  500.       
  501.     default: err(gtolger);
  502.   }
  503.   avma=av;
  504.   return y;
  505. }
  506.  
  507.  
  508. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  509. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  510. /*                                                                 */
  511. /*                    COMPARAISON  A ZERO                          */
  512. /*                                                                 */
  513. /*        x pointe sur un GEN;la fonction renvoie 1 si le          */
  514. /*                   GEN est nul,0 sinon .                         */
  515. /*                                                                 */
  516. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  517. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  518.  
  519.  
  520. gcmp0(x)
  521.      
  522.      GEN     x;
  523.      
  524. {
  525.   long    typy=typ(x),l,i;
  526.   
  527.   switch(typy)
  528.   {
  529.     case 1 :
  530.     case 2 :
  531.     case 10:
  532.     case 11: return !signe(x);
  533.       
  534.     case 3 :
  535.     case 9 : return gcmp0(x[2]);
  536.       
  537.     case 4 :
  538.     case 5 : return !signe(x[1]);
  539.     case 13:
  540.     case 14: return gcmp0(x[1]);
  541.       
  542.     case 6 : return gcmp0(x[1])&&gcmp0(x[2]);
  543.       
  544.     case 8 : return gcmp0(x[2])&&gcmp0(x[3]);
  545.       
  546.     case 7 : return !signe(x[4]);
  547.       
  548.     case 15:
  549.     case 16: 
  550.     case 17:
  551.     case 18:
  552.     case 19: l=lg(x);i=1;
  553.       while ((i<l)&&gcmp0(x[i])) i++;
  554.       return i==l;
  555.       
  556.     default: err(gcmper);
  557.   }
  558. }
  559.  
  560.  
  561. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  562. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  563. /*                                                                 */
  564. /*            COMPARAISONS  A     1  et -1                         */
  565. /*                                                                 */
  566. /*          x pointe sur un GEN;la fonction renvoie 1 si le        */
  567. /*          GEN est l'entier 1 (resp.-1),0 sinon .                 */
  568. /*                                                                 */
  569. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  570. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  571.  
  572.  
  573. gcmp1(x)
  574.      
  575.      GEN     x;
  576.      
  577. {
  578.   
  579.   long    typy=typ(x),l,y;
  580.   GEN     p1;
  581.   
  582.   switch(typy)
  583.   {
  584.     case 1 : return (lgef(x)==3)&&(signe(x)==1)&&(mant(x,1)==1);
  585.     case 2 : l=avma;p1=subrs(x,1);
  586.       y=(!signe(p1));avma=l;
  587.       return y;
  588.     case 3 :
  589.     case 9 : return gcmp1(x[2]);
  590.     case 4 : return gcmp1(x[1])&&gcmp1(x[2]);
  591.     case 5 : return !(cmpii(x[1],x[2]));
  592.     case 6 : return gcmp1(x[1])&&gcmp0(x[2]);
  593.     case 8 : return gcmp1(x[2])&&gcmp0(x[3]);
  594.     case 7 : return !valp(x)&&gcmp1(x[4]);
  595.     case 10 : return (lgef(x)==3)&&gcmp1(x[2]);
  596.     default: return 0;
  597.   }
  598. }
  599.  
  600.  
  601. gcmp_1(x)
  602.      
  603.      GEN     x;
  604.      
  605. {
  606.   
  607.   long    typy=typ(x),l,y;
  608.   GEN     p1;
  609.   
  610.   switch(typy)
  611.   {
  612.     case 1 : return (lgef(x)==3)&&(signe(x)== -1)&&(mant(x,1)==1);
  613.     case 2 : l=avma;p1=addsr(1,x);
  614.       y=(!signe(p1));avma=l;
  615.       return y;
  616.     case 3 : l=avma;p1=addsi(1,x[2]);
  617.       y=cmpii(p1,x[1]);avma=l;
  618.       return !y;
  619.     case 4 :
  620.     case 5 : l=avma;p1=addii(x[1],x[2]);
  621.       y=signe(p1);avma=l;
  622.       return !y;
  623.     case 6 : return gcmp_1(x[1])&&gcmp0(x[2]);
  624.     case 8 : return gcmp_1(x[2])&&gcmp0(x[3]);
  625.     case 7 : l=avma;p1=addsi(1,x[4]);
  626.       y=cmpii(p1,x[3]);avma=l;
  627.       return !y;
  628.     case 9 : l=avma;p1=gaddsg(1,x[2]);
  629.       y=(!gegal(p1,x[1]))&&(signe(p1));avma=l;
  630.       return !y;
  631.     case 10 : return (lgef(x)==3)&&gcmp_1(x[2]);
  632.     default: return 0;
  633.     }
  634. }
  635.  
  636. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  637. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  638. /*                                                                 */
  639. /*                       COMPARAISON SIGNEE                        */
  640. /*                                                                 */
  641. /*            rend le signe de x-y si ceci a un sens,              */
  642. /*            sinon erreur.                                        */
  643. /*                                                                 */
  644. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  645. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  646.  
  647.  
  648. long    gcmp(x,y)
  649.      
  650.      GEN     x,y;
  651.      
  652. {
  653.   long  tx,ty,f,avmacourant;
  654.   
  655.   tx=typ(x);ty=typ(y);
  656.   if ((tx>5)||(tx==3)||(ty>5)||(ty==3)) err(gcmper);
  657.   if((tx<=2)&&(ty<=2)) return mpcmp(x,y);
  658.   else
  659.   {
  660.     avmacourant=avma;f=gsigne(gsub(x,y));
  661.     avma=avmacourant;return f;
  662.   }
  663. }
  664.  
  665. long lexcmp(x,y)
  666.      GEN x,y;
  667. {
  668.   long tx=typ(x),ty=typ(y),lx,ly,fl,s,i;
  669.   GEN p1;
  670.  
  671.   if(tx>ty) {p1=x;x=y;y=p1;fl=tx;tx=ty;ty=fl;s= -1;} else s=1;
  672.   ly=lg(y);
  673.   if(tx<17)
  674.     {
  675.       if(ty<17) return s*gcmp(x,y);
  676.       if(ly==1) return s;
  677.       fl=lexcmp(x,y[1]);if(fl) return s*fl;
  678.       return (ly>2)?-s:0;
  679.     }
  680.   lx=lg(x);
  681.   if(ly==1) return (lx==1)?0:s;
  682.   if(lx==1) return -s;
  683.   if((ty==19)&&(tx<19))
  684.     {
  685.       fl=lexcmp(x,y[1]);if(fl) return s*fl;
  686.       return (ly>2)?-s:0;
  687.     }
  688.   if(lx>ly) {p1=x;x=y;y=p1;fl=lx;lx=ly;ly=fl;s= -s;}
  689.   fl=0;for(i=1;(i<lx)&&(!fl);i++) fl=lexcmp(x[i],y[i]);
  690.   if(fl) return s*fl;
  691.   return (ly>lx)?-s:0;
  692. }
  693.  
  694. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  695. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  696. /*                                                                 */
  697. /*                            EGALITE                              */
  698. /*                                                                 */
  699. /*            Renvoie 1 si x=y, 0 sinon                            */
  700. /*                                                                 */
  701. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  702. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  703.  
  704.  
  705. long gegal(x,y)
  706.      
  707.      GEN     x,y;
  708.      
  709. {
  710.   long    avmacourant,f,i,tx=typ(x),ty=typ(y),lx,masq=0xff00ffff;
  711.   GEN     p1;
  712.   
  713.   if(tx!=ty)
  714.   {
  715.     avmacourant=avma;p1=gsub(x,y);f=gcmp0(p1);avma=avmacourant;
  716.   }
  717.   else
  718.   {
  719.     switch(tx)
  720.     {
  721.       case 1: if((x[1]&masq)!=(y[1]&masq)) return 0;
  722.         i=2;lx=lgef(x);while((i<lx)&&(x[i]==y[i])) i++;return i==lx;
  723.       case 10: if(x[1]!=y[1]) return 0;
  724.         i=2;lx=lgef(x);while((i<lx)&&(gegal(x[i],y[i]))) i++;return i==lx;
  725.       case 3:
  726.       case 6:
  727.       case 9: return gegal(x[1],y[1])&&gegal(x[2],y[2]);
  728.       case 8: return gegal(x[1],y[1])&&gegal(x[2],y[2])&&gegal(x[3],y[3]);
  729.       case 4:
  730.       case 5:
  731.       case 13:
  732.       case 14: avmacourant=avma;f=gegal(gmul(x[1],y[2]),gmul(x[2],y[1]));
  733.     avma=avmacourant;return f;
  734.       default: avmacourant=avma;p1=gsub(x,y);f=gcmp0(p1);avma=avmacourant;
  735.     }
  736.   }
  737.   return f;
  738. }
  739.  
  740. long polegal(x,y)
  741.      
  742.      GEN    x,y;
  743.      
  744.      /* a usage interne donc pas de verification de types */
  745.      
  746. {
  747.   long i;
  748.   
  749.   if (x[1] != y[1]) return 0;
  750.   for(i = lgef(x) - 1; i > 1; i--) if (!gegal(x[i],y[i])) return 0;
  751.   return 1;
  752. }
  753.  
  754. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  755. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  756. /*                                                                 */
  757. /*                          VALUATION                              */
  758. /*                                                                 */
  759. /*   renvoie le plus grand exposant de p divisant x quand cela a   */
  760. /*   un sens (erreur pour les types reels, integermod et polymod   */
  761. /*   si le module n'est pas divisible par p, et q-adiques pour     */
  762. /*   q!=p. p doit etre de type entier ou polynome.                 */
  763. /*                                                                 */
  764. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  765. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  766.  
  767. long ggval(x,p)
  768.      
  769.      GEN x,p;
  770.      
  771. {
  772.   long  tx=typ(x),tp=typ(p),lx,i,j,vx,v,av,val,limite=(avma + bot) / 2, tetpil;
  773.   GEN  p1,p2,reste;
  774.   
  775.   if(isexactzero(x)) err(gvaler2);
  776.   if((tx<9)&&(tp==10)) return 0;
  777.   switch(tx)
  778.     {
  779.     case 1: if(tp!=1) err(gvaler4);
  780.       val=0;av=avma;
  781.       while (1)
  782.     {
  783.       p1=dvmdii(x,p,&reste);
  784.       if (signe(reste)) break;
  785.       val++;x=p1;
  786.           if (avma < limite)
  787.           {
  788.             tetpil = avma;
  789.             x = gerepile(av, tetpil, gcopy(x));
  790.           }
  791.     }
  792.       avma=av;break;
  793.     case 3: p1=cgeti(lgef(x[1]));
  794.       if((tp!=1)||(!mpdivis(x[1],p,p1))) err(gvaler4);
  795.       p2=cgeti(lgef(x[2]));
  796.       if(mpdivis(x[2],p,p2))
  797.     {val=1;while(mpdivis(p1,p,p1)&&mpdivis(p2,p,p2)) val++;}
  798.       else val=0;
  799.       break;
  800.     case 7: if((tp!=1)||(!gegal(p,x[2]))) err(gvaler4);
  801.       val=valp(x);break;
  802.     case 9: if(tp==1) val=ggval(x[2],p);
  803.       else
  804.     {
  805.       if((tp!=10)||(!poldivis(x[1],p,p1))) err(gvaler4);
  806.       if(poldivis(x[2],p,p2))
  807.         {val=1;while(poldivis(p1,p,p1)&&poldivis(p2,p,p2)) val++;}
  808.       else val=0;
  809.     }
  810.       break;
  811.     case 10: i=2;while (isexactzero(x[i])) i++;
  812.       if(tp==10)
  813.     {
  814.       v=varn(p);vx=varn(x);if(vx>v) return 0;
  815.       if(vx==v)
  816.         {
  817.           if(ismonome(p)) val=i-2;
  818.           else
  819.         {
  820.           val=0;av=avma;
  821.           while (1)
  822.             {
  823.               p1=poldivres(x,p,&reste);
  824.               if (signe(reste)) break;
  825.               val++;x=p1;
  826.                       if (avma < limite)
  827.                       {
  828.                         tetpil = avma;
  829.                         x = gerepile(av, tetpil, gcopy(x));
  830.                       }
  831.             }
  832.           avma=av;
  833.         }
  834.         }
  835.       else
  836.         {
  837.           val=ggval(x[i],p);
  838.           for(j=i+1;j<lgef(x);j++)
  839.         if(!gcmp0(x[j])) val=min(val,ggval(x[j],p));
  840.         }
  841.     }
  842.       else 
  843.     {
  844.       if(tp!=1) err(gvaler4);
  845.       val=ggval(x[i],p);
  846.       for(j=i+1;j<lgef(x);j++)
  847.         if(!gcmp0(x[j])) val=min(val,ggval(x[j],p));
  848.     }
  849.       break;
  850.       
  851.     case 11: if((tp!=10)&&(tp!=11)&&(tp!=1)) err(gvaler4);
  852.       v=gvar(p);vx=varn(x);if(vx>v) return 0;
  853.       if(vx==v) return (long)(valp(x)/ggval(p,polx[v]));
  854.       val=ggval(x[2],p);
  855.       for(j=3;j<lg(x);j++)
  856.         if(!gcmp0(x[j])) val=min(val,ggval(x[j],p));
  857.       break;
  858.     case 4:
  859.     case 5:
  860.     case 13:
  861.     case 14: val=ggval(x[1],p)-ggval(x[2],p);break;
  862.       
  863.     case 2:
  864.     case 15:
  865.     case 16: err(gvaler4);
  866.     case 6:
  867.     case 8:
  868.     case 17:
  869.     case 18:
  870.     case 19: val=VERYBIGINT;lx=lg(x);
  871.       for(j=1;j<lx;j++) val=min(val,ggval(x[j],p));
  872.       break;
  873.     default: err(gvaler4);
  874.   }
  875.   return val;
  876. }
  877.  
  878. pvaluation(x,p,py)
  879.      
  880.      GEN     x,p,*py;
  881.      
  882. {
  883.   long l,tetpil,val,limite = (bot + avma) / 2;
  884.   GEN  p1,reste;
  885.   
  886.   val=0;*py=cgeti(lgef(x));l=avma;
  887.   while (1)
  888.   {
  889.     p1=dvmdii(x,p,&reste);
  890.     if (signe(reste)) break;
  891.     val++;x=p1;
  892.     if (avma < limite)
  893.       {
  894.         tetpil = avma;
  895.         x = gerepile(l, tetpil, gcopy(x));
  896.       }
  897.   }
  898.   affii(x,*py);avma=l;
  899.   return val;
  900. }
  901.  
  902.  
  903. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  904. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  905. /*                                                                 */
  906. /*                            NEGATION                             */
  907. /*                                                                 */
  908. /*            Cree-x,ou x pointe sur un GEN .                      */
  909. /*                                                                 */
  910. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  911. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  912.  
  913.  
  914. GEN     gneg(x)
  915.      
  916.      GEN     x;
  917.      
  918. {
  919.   long  tx,lx,i,sx;
  920.   GEN   y;
  921.   
  922.   if (gcmp0(x))
  923.   {
  924.     y=gcopy(x);
  925.   }
  926.   else
  927.   {
  928.     tx=typ(x);lx=lg(x);
  929.     
  930.     switch(tx)
  931.     {
  932.       case 1 :
  933.       case 2 : y=gcopy(x);
  934.         sx=signe(x);
  935.         sx= -sx;
  936.         setsigne(y,sx);
  937.         break;
  938.         
  939.       case 3 : y=cgetg(3,tx);y[1]=copyifstack(x[1]);
  940.         if(gcmp0(x[2])) y[2]=zero;
  941.         else
  942.           y[2]=lsubii(y[1],x[2]);
  943.         break;
  944.         
  945.       case 9: y=cgetg(3,tx);y[1]=copyifstack(x[1]);
  946.         y[2]=lneg(x[2]);
  947.         break;
  948.         
  949.       case 4 :
  950.       case 5 :
  951.       case 13:
  952.       case 14: y=cgetg(3,tx);y[1]=lneg(x[1]);y[2]=lcopy(x[2]);
  953.         break;
  954.         
  955.       case 7 : y=cgetp(x);
  956.         if(gcmp0(x[4])) {affsi(0,y[4]);y[1]=x[1];}
  957.         else
  958.           {
  959.             subiiz(y[3],x[4],y[4]);
  960.             setvalp(y,valp(x));
  961.           }
  962.         break;
  963.  
  964.       case 8 : y=cgetg(lx,tx);
  965.         for (i=2;i<lx;i++) y[i]=lneg(x[i]);
  966.         y[1]=copyifstack(x[1]);
  967.         break;
  968.         
  969.       case 6 :
  970.       case 17:
  971.       case 18:
  972.       case 19: y=cgetg(lx,tx);
  973.         for (i=1;i<lx;i++) y[i]=lneg(x[i]);
  974.         break;
  975.  
  976.       case 15:
  977.       case 16: err(gneger);
  978.         
  979.       case 10: lx=lgef(x);
  980.         y=cgetg(lx,tx);
  981.         for (i=2;i<lx;i++) y[i]=lneg(x[i]);
  982.         y[1]=x[1];
  983.         break;
  984.         
  985.       case 11: y=cgetg(lx,tx);
  986.         for (i=2;i<lx;i++) y[i]=lneg(x[i]);
  987.         y[1]=x[1];
  988.     }
  989.   }
  990.   return y;
  991. }
  992.  
  993.  
  994. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  995. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  996. /*                                                                 */
  997. /*                        VALEUR ABSOLUE                           */
  998. /*                                                                 */
  999. /*            Cree un GEN pointant sur la valeur absolue de x      */
  1000. /*            a condition que x soit un entier,ou un reel,         */
  1001. /*            ou une fraction,ou un complexe;sinon erreur .        */
  1002. /*                                                                 */
  1003. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1004. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1005.  
  1006. GEN     gabs(x,prec)
  1007.      
  1008.      GEN     x;
  1009.      
  1010. {
  1011.   long  tx=typ(x),lx=lg(x),i,l,tetpil;
  1012.   GEN   y,p1;
  1013.   
  1014.   switch(tx)
  1015.   {
  1016.     case 1 :
  1017.     case 2 : y=mpabs(x);break;
  1018.       
  1019.     case 4 :
  1020.     case 5 : y=cgetg(lx,tx);
  1021.       y[1]=lmpabs(x[1]);
  1022.       y[2]=lmpabs(x[2]);
  1023.       break;
  1024.       
  1025.     case 6 : l=avma;p1=gnorm(x);tetpil=avma;
  1026.       y=gsqrt(p1,prec);y=gerepile(l,tetpil,y);
  1027.       break;
  1028.       
  1029.     case 8 : y=transc(gabs,x,prec);break;
  1030.       
  1031.     case 17:
  1032.     case 18:
  1033.     case 19: y=cgetg(lx,tx);
  1034.       for(i=1;i<lx;i++)
  1035.         y[i]=labs(x[i]);
  1036.       break;
  1037.       
  1038.     default: err(gabser);
  1039.   }
  1040.   return y;
  1041. }
  1042.  
  1043. GEN     gmax(x,y)
  1044.      GEN   x,y;
  1045.      
  1046. {
  1047.   return (gcmp(x,y)>=0) ? gcopy(x) : gcopy(y);
  1048. }
  1049.  
  1050. GEN     gmin(x,y)
  1051.      GEN   x,y;
  1052.      
  1053. {
  1054.   return (gcmp(x,y)<=0) ? gcopy(x) : gcopy(y);
  1055. }
  1056.  
  1057.  
  1058. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1059. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1060. /*                                                                 */
  1061. /*                  AFFECTATION  S--> GEN                          */
  1062. /*                                                                 */
  1063. /*            Etant donnes un long et un GEN,affecte le long       */
  1064. /*            dans le GEN,permettant une initialisation .          */
  1065. /*                                                                 */
  1066. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1067. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1068.  
  1069.  
  1070. void gaffsg(s,x)
  1071.      
  1072.      long    s;
  1073.      GEN     x;
  1074.      
  1075. {
  1076.   long  tx,l,i,v,val;
  1077.   GEN   reste,p1;
  1078.   
  1079.   tx=typ(x);
  1080.   
  1081.   switch(tx)
  1082.   {
  1083.     case 1 : affsi(s,x);break;
  1084.     case 2 : affsr(s,x);break;
  1085.     case 3 : modsiz(s,x[1],x[2]);break;
  1086.     case 4 :
  1087.     case 5 : affsi(s,x[1]);
  1088.       affsi(1,x[2]);
  1089.       break;
  1090.     case 6 : gaffsg(s,x[1]);
  1091.       gaffsg(0,x[2]);
  1092.       break;
  1093.     case 7 : if (!s) {setvalp(x,defaultpadicprecision);x[4]=zero;}
  1094.     else
  1095.     {
  1096.       val=0;
  1097.       l=avma;
  1098.       while (1)
  1099.       {
  1100.         p1=dvmdsi(s,x[2],&reste);
  1101.         if (signe(reste)) break;
  1102.         val++;s=itos(p1);
  1103.       }
  1104.       avma=l;
  1105.       setvalp(x,val);
  1106.       if (s>0) affsi(s,x[4]);
  1107.       else addsiz(s,x[3],x[4]);
  1108.     }
  1109.     break;
  1110.       
  1111.     case 8 : gaffsg(s,x[2]);
  1112.       gaffsg(0,x[3]);
  1113.       break;
  1114.       
  1115.     case 9 : gaffsg(s,x[2]);break;
  1116.       
  1117.     case 10: v=varn(x);
  1118.       if (!s) x[1]=2;
  1119.       else
  1120.       {
  1121.         x[1]=0x1000003;gaffsg(s,x[2]);
  1122.       }
  1123.       setvarn(x,v);break;
  1124.       
  1125.     case 11: v=varn(x);gaffsg(s,x[2]);l=lg(x);
  1126.       if (!s) x[1]=0x7ffe +l;
  1127.       else
  1128.       {
  1129.         x[1]=0x1008000;
  1130.         for (i=3;i<l;gaffsg(0,x[i]),i++);
  1131.       }
  1132.       setvarn(x,v);break;
  1133.       
  1134.     case 13:
  1135.     case 14: gaffsg(s,x[1]);gaffsg(1,x[2]);
  1136.       break;
  1137.       
  1138.     case 17:
  1139.     case 18: if (lg(x)!=2) err(gaffsger1);
  1140.     else gaffsg(s,x[1]);break;
  1141.       
  1142.     case 19: if (lg(x)!=2) err(gaffsger2);
  1143.     else gaffsg(s,x[1]);break;
  1144.       
  1145.     default: err(gaffer1);
  1146.   }
  1147. }
  1148.  
  1149.  
  1150. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1151. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1152. /*                                                                 */
  1153. /*                    AFFECTATION GENERALE                         */
  1154. /*                                                                 */
  1155. /*            Etant donnes deux GEN x et y,affecte le contenu      */
  1156. /*            de x dans y,si cela est possible .                   */
  1157. /*                                                                 */
  1158. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1159. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1160.  
  1161.  
  1162. void gaffect(x,y)
  1163.      
  1164.      GEN     x,y;
  1165.      
  1166. {
  1167.   long  i,j,k,l,v,vy,val,lx,ly,tx,ty,d,avmacourant;
  1168.   long  num,den;
  1169.   GEN   p1;
  1170.   
  1171.   tx=typ(x);ty=typ(y);lx=lg(x);ly=lg(y);
  1172.   
  1173.   
  1174.   if (tx<10)
  1175.   {
  1176.     if (ty>=10)
  1177.     {
  1178.       switch(ty)
  1179.       {
  1180.         case 10: v=varn(y);gaffect(x,y[2]);
  1181.           for(i=3;i<ly;gaffsg(0,y[i]),i++);
  1182.           if (gcmp0(x)) y[1]=2;
  1183.           else y[1]=0x1000003;
  1184.           setvarn(y,v);break;
  1185.           
  1186.         case 11: v=varn(y);gaffect(x,y[2]);
  1187.           if (gcmp0(x)) y[1]=0x7ffe +ly;
  1188.           else
  1189.           {
  1190.             y[1]=0x1008000;
  1191.             for (i=3;i<ly;gaffsg(0,y[i]),i++);
  1192.           }
  1193.           setvarn(y,v);break;
  1194.           
  1195.         case 13:
  1196.         case 14: gaffect(x,y[1]);gaffsg(1,y[2]); break;
  1197.         case 17:
  1198.         case 18:
  1199.         case 19: err(gaffer2);
  1200.           
  1201.         default: err(gaffer1);
  1202.         }
  1203.     }
  1204.     else
  1205.     {
  1206.       switch(tx)
  1207.       {
  1208.         case 1 :
  1209.           switch(ty)
  1210.           {
  1211.             case 1 : affii(x,y);break;
  1212.             case 2 : affir(x,y);break;
  1213.             case 3 : modiiz(x,y[1],y[2]);break;
  1214.             case 4 :
  1215.             case 5 : affii(x,y[1]);affsi(1,y[2]); break;
  1216.             case 6 : gaffect(x,y[1]);gaffsg(0,y[2]); break;
  1217.             case 7 : if (!signe(x)) {setvalp(y,defaultpadicprecision);y[4]=zero;}
  1218.               else
  1219.               {
  1220.                 l=avma;
  1221.                 val=pvaluation(x,y[2],&p1);
  1222.                 setvalp(y,val);modiiz(p1,y[3],y[4]);
  1223.                 avma=l;
  1224.               }
  1225.               break;
  1226.             case 8 : gaffect(x,y[2]);gaffsg(0,y[3]); break;
  1227.         case 9 : gaffect(x,y[2]);break;
  1228.             default: err(gaffer1);
  1229.             }
  1230.             break;
  1231.           
  1232.         case 2 :
  1233.           switch(ty)
  1234.           {
  1235.             case 1 : err(gaffer3);
  1236.             case 2 : affrr(x,y);break;
  1237.             case 3 :
  1238.             case 4 :
  1239.             case 5 : err(gaffer3);
  1240.             case 6 : gaffect(x,y[1]);gaffsg(0,y[2]); break;
  1241.             case 7 :
  1242.             case 8 : err(gaffer3);
  1243.         case 9 : gaffect(x,y[2]);break;
  1244.             default: err(gaffer1);
  1245.           }
  1246.           break;
  1247.           
  1248.         case 3 :
  1249.           switch(ty)
  1250.           {
  1251.             case 1 :
  1252.             case 2 : err(gaffer4);
  1253.             case 3 : if (!divise(x[1],y[1])) err(gaffer5);
  1254.               modiiz(x[2],y[1],y[2]); break;
  1255.             case 4 :
  1256.             case 5 :
  1257.             case 6 :
  1258.             case 8 : err(gaffer4);
  1259.             case 7 : err(gaffer6);
  1260.         case 9 : gaffect(x,y[2]);break;
  1261.             default: err(gaffer1);
  1262.           }
  1263.           break;
  1264.         case 4 :
  1265.           switch(ty)
  1266.           {
  1267.             case 1 : i=mpdivis(x[1],x[2],y);
  1268.               if (!i) err(gaffer7);
  1269.               break;
  1270.             case 2 : diviiz(x[1],x[2],y);break;
  1271.               
  1272.             case 3 : avmacourant=avma;
  1273.               p1=mpinvmod(x[2],y[1]);
  1274.               modiiz(mulii(x[1],p1),y[1],y[2]);
  1275.               avma=avmacourant;
  1276.               break;
  1277.             case 4 :
  1278.             case 5 : for (i=1;i<=2;affii(x[i],y[i]),i++);break;
  1279.             case 6 : gaffect(x,y[1]);gaffsg(0,y[2]); break;
  1280.             case 7 : if(!signe(x[1]))
  1281.               {
  1282.                 setvalp(y,defaultpadicprecision);y[4]=zero;
  1283.               }
  1284.               else
  1285.               {
  1286.                 l=avma;num=x[1];den=x[2];
  1287.                 if(!(val=pvaluation(num,y[2],&num)))
  1288.                   val= -pvaluation(den,y[2],&den);
  1289.                 p1=mulii(num,mpinvmod(den,y[3]));
  1290.                 modiiz(p1,y[3],y[4]);avma=l;
  1291.                 setvalp(y,val);
  1292.               }
  1293.               break;
  1294.             case 8 : gaffect(x,y[2]);gaffsg(0,y[3]); break;
  1295.         case 9 : gaffect(x,y[2]);break;
  1296.             default: err(gaffer1);
  1297.           }
  1298.           break;
  1299.         case 5 :
  1300.           switch(ty)
  1301.           {
  1302.             case 1 : i=mpdivis(x[1],x[2],y);
  1303.               if (!i) err(gaffer7);
  1304.               break;
  1305.             case 2 : diviiz(x[1],x[2],y);break;
  1306.               
  1307.             case 3 : avmacourant=avma;gaffect(gred(x),y);
  1308.               avma=avmacourant;
  1309.               break;
  1310.             case 4 : gredz(x,y);break;
  1311.               
  1312.             case 5 : for (i=1;i<=2;affii(x[i],y[i]),i++);break;
  1313.             case 6 : gaffect(x,y[1]);gaffsg(0,y[2]); break;
  1314.             case 7 : if(!signe(x[1]))
  1315.               {
  1316.                 setvalp(y,defaultpadicprecision);y[4]=zero;
  1317.               }
  1318.               else
  1319.               {
  1320.                 l=avma;num=x[1];den=x[2];
  1321.                 val=pvaluation(num,y[2],&num)-pvaluation(den,y[2],&den);
  1322.                 p1=mulii(num,mpinvmod(den,y[3]));
  1323.                 modiiz(p1,y[3],y[4]);avma=l;
  1324.                 setvalp(y,val);
  1325.               }
  1326.               break;
  1327.             case 8 : gaffect(x,y[2]);gaffsg(0,y[3]); break;
  1328.         case 9 : gaffect(x,y[2]);break;
  1329.             default: err(gaffer1);
  1330.             }
  1331.             break;
  1332.             
  1333.         case 6 :
  1334.           switch(ty)
  1335.           {
  1336.             case 1 :
  1337.             case 2 :
  1338.             case 3 :
  1339.             case 4 :
  1340.             case 5 :
  1341.             case 7 :
  1342.             case 8 : if (!gcmp0(x[2])) err(gaffer8);
  1343.               gaffect(x[1],y);
  1344.               break;
  1345.               
  1346.             case 6 : for (i=1;i<=2;gaffect(x[i],y[i]),i++);break;
  1347.         case 9 : gaffect(x,y[2]);break;
  1348.             default: err(gaffer1);
  1349.           }
  1350.           break;
  1351.           
  1352.         case 7 :
  1353.           switch(ty)
  1354.           {
  1355.             case 1 :
  1356.             case 2 : err(gaffer10);
  1357.             case 3 : if(valp(x)<0) err(gaffer11);
  1358.               l=avma;
  1359.               val=pvaluation(y[1],x[2],&p1);
  1360.               if(!gcmp1(p1)) err(gaffer11);
  1361.               if(val>(signe(x[4])?(precp(x)+valp(x)):valp(x)+1)) err(gaffer11);
  1362.               modiiz(x[4],y[1],y[2]);
  1363.               avma=l;break;
  1364.             case 4 :
  1365.             case 5 :
  1366.             case 6 :
  1367.             case 8 : err(gaffer10);
  1368.               
  1369.             case 7 : if(cmpii(x[2],y[2])) err(gaffer12);
  1370.               modiiz(x[4],y[3],y[4]);
  1371.               setvalp(y,valp (x));
  1372.               break;
  1373.               
  1374.         case 9 : gaffect(x,y[2]);break;
  1375.             default: err(gaffer1);
  1376.           }
  1377.           break;
  1378.           
  1379.         case 8 :
  1380.           switch(ty)
  1381.           {
  1382.             case 1 :
  1383.             case 3 :
  1384.             case 4 :
  1385.             case 5 :
  1386.             case 7 : if(!gcmp0(x[3])) err(gaffer9);
  1387.               gaffect(x[2],y);
  1388.               break;
  1389.             case 2 : l=avma;p1=co8(x,ly);gaffect(p1,y);
  1390.               avma=l;break;
  1391.             case 6 : ly=precision(y);
  1392.               if(ly)
  1393.               {
  1394.                 l=avma;p1=co8(x,ly);gaffect(p1,y);
  1395.                 avma=l;
  1396.               }
  1397.               else
  1398.               {
  1399.                 if(!gcmp0(x[3])) err(gaffer9);
  1400.                 gaffect(x[2],y);
  1401.               }
  1402.               break;
  1403.             case 8 : if(!gegal(x[1],y[1])) err(gaffer21);
  1404.               for (i=2;i<=3;gaffect(x[i],y[i]),i++);
  1405.               break;
  1406.         case 9 : gaffect(x,y[2]);break;
  1407.             default: err(gaffer1);
  1408.           }
  1409.           break;
  1410.     case 9 : if(ty!=9) err(gaffer17);
  1411.       if (gdivise(x[1],y[1]))
  1412.         gmodz(x[2],y[1],y[2]);else err(gaffer18);
  1413.       break;
  1414.         default: err(gaffer1);
  1415.     }
  1416.     }
  1417.   }
  1418.   else
  1419.   {
  1420.     if (ty<9) err(gaffer13);
  1421.     switch(tx)
  1422.     {
  1423.       case 10: v=varn(x);
  1424.         switch(ty)
  1425.         {
  1426.           case 10: if((vy=varn(y))>v) err(gaffer13);
  1427.             if(vy==v)
  1428.             {
  1429.               d=lgef(x);
  1430.               if (d>ly) err(gaffer14);
  1431.               y[1]=x[1];
  1432.               for (i=2;i<d;gaffect(x[i],y[i]),i++);
  1433.             }
  1434.             else
  1435.             {
  1436.               gaffect(x,y[2]);
  1437.               for(i=3;i<ly;gaffsg(0,y[i]),i++);
  1438.               if (signe(x)) y[1]=0x1000003;
  1439.               else y[1]=2;
  1440.               setvarn(y,vy);
  1441.             }
  1442.             break;
  1443.             
  1444.           case 11: if((vy=varn(y))>v) err(gaffer13);
  1445.             if (!signe(x)) gaffsg(0,y);
  1446.             else
  1447.             {
  1448.               if(vy==v)
  1449.               {
  1450.                 i=gval(x,v);setvalp(y,i);setsigne(y,1);
  1451.                 k=lgef(x)-i;
  1452.                 if(k>ly) k=ly;
  1453.                 for (j=2;j<k;j++)
  1454.                   gaffect(x[i+j],y[j]);
  1455.                 for (j=k;j<ly;j++)
  1456.                   gaffsg(0,y[j]);
  1457.               }
  1458.               else
  1459.               {
  1460.                 gaffect(x,y[2]);
  1461.                 if (!signe(x)) y[1]=0x7ffe +ly;
  1462.                 else
  1463.                 {
  1464.                   y[1]=0x1008000;
  1465.                   for (i=3;i<ly;gaffsg(0,y[i]),i++);
  1466.                 }
  1467.                 setvarn(y,vy);
  1468.               }
  1469.             }
  1470.             break;
  1471.             
  1472.           case 9: gmodz(x,y[1],y[2]);break;
  1473.           case 13:
  1474.           case 14: gaffect(x,y[1]);gaffsg(1,y[2]); break;
  1475.           case 17:
  1476.           case 18:
  1477.           case 19: err(gaffer15);
  1478.             
  1479.           default: err(gaffer1);
  1480.         }
  1481.         break;
  1482.         
  1483.       case 11: if (ty!=11) err(gaffer16);
  1484.         v=varn(x);if((vy=varn(y))>v) err(gaffer13);
  1485.         if(vy==v)
  1486.         {
  1487.           y[1]=x[1];
  1488.           if(!gcmp0(x))
  1489.           {
  1490.             k=lx;if (k>ly) k=ly;
  1491.             for (i=2;i<k;i++) gaffect(x[i],y[i]);
  1492.             for (i=k;i<ly;i++) gaffsg(0,y[i]);
  1493.           }
  1494.         }
  1495.         else
  1496.         {
  1497.           gaffect(x,y[2]);
  1498.           if (!signe(x)) y[1]=0x7ffe +ly;
  1499.           else
  1500.           {
  1501.             y[1]=0x1008000;
  1502.             for (i=3;i<ly;gaffsg(0,y[i]),i++);
  1503.           }
  1504.           setvarn(y,vy);
  1505.         }
  1506.         break;
  1507.         
  1508.       case 13:
  1509.       case 14: switch(ty)
  1510.       {
  1511.         case 10:
  1512.         case 17:
  1513.         case 18:
  1514.         case 19: err(gaffer19);
  1515.         case 11: gdivz(x[1],x[2],y);break;
  1516.         case 9: avmacourant=avma;p1=ginvmod(x[2],y[1]);
  1517.           gmod(gmul(x[1],p1),y[1],y[2]);
  1518.           avma=avmacourant;
  1519.           break;
  1520.           
  1521.         case 13:
  1522.           if (tx==14) gredz(x,y);
  1523.           else
  1524.             for (i=1;i<=2;gaffect(x[i],y[i]),i++);
  1525.           break;
  1526.           
  1527.         case 14: for (i=1;i<=2;gaffect(x[i],y[i]),i++); break;
  1528.           
  1529.         default: err(gaffer1);
  1530.       }
  1531.       break;
  1532.  
  1533.       case 15:
  1534.       case 16:
  1535.       case 17:
  1536.       case 18:
  1537.       case 19: if (ty!=tx) err(gaffer20);
  1538.         if (ly!=lx) err(gaffer20);
  1539.         for (i=1;i<lx;gaffect(x[i],y[i]),i++);
  1540.         break;
  1541.         
  1542.       default: err(gaffer1);
  1543.     }
  1544.   }
  1545. }
  1546.  
  1547. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1548. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1549. /*                                                                 */
  1550. /*                 CONVERSION QUAD-->REEL,COMPLEXE                 */
  1551. /*                          OU P-ADIQUE                            */
  1552. /*                                                                 */
  1553. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1554. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1555.  
  1556. GEN     co8(x,l)
  1557.      
  1558.      GEN     x;
  1559.      long    l;
  1560.      
  1561. {
  1562.   long  av,tetpil;
  1563.   GEN   y,p1,p2,p3;
  1564.   
  1565.   p3=(GEN)x[1];
  1566.   av=avma;p1=gmul2n(p3[3],-1);
  1567.   p2=gsqrt(gsub(gmul(p1,p1),p3[2]),l);
  1568.   p1=gmul(x[3],gsub(p2,p1));
  1569.   tetpil=avma;y=gerepile(av,tetpil,gadd(p1,x[2]));
  1570.   return y;
  1571. }
  1572.  
  1573. GEN   cvtop(x,p,l)
  1574.      
  1575.      GEN  x,p;
  1576.      long l;
  1577.      
  1578. {
  1579.   GEN z,p1,p2,p3;
  1580.   long av,tetpil,n,tx=typ(x);
  1581.   
  1582.   if(typ(p)!=1) err(cvtoper1);
  1583.   if(gcmp0(x)) z=ggrando(p,l);
  1584.   else
  1585.   {
  1586.     switch(tx)
  1587.     {
  1588.       case 1: z=gadd(x,ggrando(p,ggval(x,p)+l));break;
  1589.       case 2: err(cvtoper2);
  1590.       case 3: n=ggval(x[1],p);
  1591.         if(n>l) n=l;
  1592.         z=gadd(x[2],ggrando(p,n));break;
  1593.       case 4:
  1594.       case 5: n=ggval(x[1],p);
  1595.         n-=ggval(x[2],p);
  1596.         z=gadd(x,ggrando(p,n+l));break;
  1597.       case 6: av=avma;p1=gsqrt(gaddgs(ggrando(p,l),-1));
  1598.         p1=gmul(p1,x[2]);tetpil=avma;
  1599.         z=gerepile(av,tetpil,gadd(p1,x[1]));
  1600.         break;
  1601.       case 7: z=gprec(x,l);break;
  1602.       case 8:
  1603.         av=avma;p1=(GEN)x[1];p3=gmul2n(p1[3],-1);
  1604.         p2=gsub(gmul(p3,p3),p1[2]);
  1605.         switch(typ(p2))
  1606.           {
  1607.           case 1: n=ggval(p2,p);
  1608.             p2=gadd(p2,ggrando(p,n+l));break;
  1609.           case 3: break;
  1610.           case 4:
  1611.           case 5: n=ggval(p2[1],p);
  1612.             n-=ggval(p2[2],p);
  1613.             p2=gadd(p2,ggrando(p,n+l));break;
  1614.           case 7: break;
  1615.           default: err(gadder14);
  1616.           }
  1617.         p2=gsqrt(p2);
  1618.         p1=gmul(x[3],gsub(p2,p3));tetpil=avma;
  1619.         z=gerepile(av,tetpil,gadd(x[2],p1));
  1620.         break;
  1621.       default: err(cvtoper2);
  1622.     }
  1623.   }
  1624.   return z;
  1625. }
  1626.  
  1627. GEN gcvtop(x,p,r)
  1628.      GEN x,p;
  1629.      long r;
  1630.  
  1631. {
  1632.   long i,tx=typ(x),lx;
  1633.   GEN y;
  1634.  
  1635.   if(tx<9) return cvtop(x,p,r);
  1636.   switch(tx)
  1637.     {
  1638.     case 10: lx=lgef(x);y=cgetg(lx,10);y[1]=x[1];
  1639.       for(i=2;i<lx;i++) y[i]=(long)gcvtop(x[i],p,r);break;
  1640.     case 11:
  1641.       if(gcmp0(x)) y=gcopy(x);
  1642.       else 
  1643.     {
  1644.       lx=lg(x);y[1]=x[1];
  1645.       for(i=2;i<lx;i++) y[i]=(long)gcvtop(x[i],p,r);
  1646.     }
  1647.       break;
  1648.     case 9:
  1649.     case 13:
  1650.     case 14:
  1651.     case 17:
  1652.     case 18:
  1653.     case 19: lx=lg(x);y=cgetg(lx,tx);
  1654.       for(i=1;i<lx;i++) y[i]=(long)gcvtop(x[i]);break;
  1655.     default: err(cvtoper2);
  1656.     }
  1657.   return y;
  1658. }
  1659.  
  1660. long   gexpo(x)
  1661.      GEN   x;
  1662.      
  1663. {
  1664.   long tx=typ(x),lx=lg(x),e,i,y,av;
  1665.   GEN  p1;
  1666.   
  1667.   switch(tx)
  1668.   {
  1669.     case 1 :if(!signe(x)) err(gexpoer2);
  1670.       return expi(x);
  1671.     case 4 :
  1672.     case 5 : if(!signe(x[1])) err(gexpoer2);
  1673.       return expi(x[1])-expi(x[2]);
  1674.     case 2 : return expo(x);
  1675.     case 6 : av=avma;p1=gnorm(x);y=(gexpo(p1)>>1); avma=av;break;
  1676.     case 8 :  if(gcmp0(x)) err(gexpoer2);
  1677.       av=avma;p1=cgetg(3,6);p1[1]=lgetr(3);
  1678.       p1[2]=lgetr(3);gaffect(x,p1);y=gexpo(p1);
  1679.       avma=av;break;
  1680.     case 10: lx=lgef(x);
  1681.     case 11: 
  1682.     case 17:
  1683.     case 18:
  1684.     case 19: y= -BIGINT;
  1685.       for(i=lontyp[tx];i<lx;i++)
  1686.       {
  1687.         e=gexpo(x[i]);if(e>y) y=e;
  1688.       }
  1689.       break;
  1690.     default: err(gexpoer1);
  1691.   }
  1692.   return y;
  1693. }
  1694.  
  1695. long gsize(x)
  1696.      GEN x;
  1697. {
  1698.   long l;
  1699.   if(gcmp0(x)) return 0;
  1700.   l=gexpo(x)*L2SL10;
  1701.   return l;
  1702. }
  1703.  
  1704. void normalize(px)
  1705.      
  1706.      GEN  *px;
  1707.      
  1708. {
  1709.   long    i,j,v,l,tetpil,e,lx,f;
  1710.   GEN     p1,x;
  1711.   
  1712.   if(typ(x= *px)!=11) err(gnormaler);
  1713.   i=2;f=1;lx=lg(x);v=varn(x);
  1714.   while (f&&(i<lx))
  1715.     {f=isexactzero(x[i]);i++;}
  1716. /*    {f=gcmp0(x[i]);i++;} */
  1717.   i--;
  1718.   if(i>2)
  1719.   {
  1720.     l=(long)(x+lx);e=valp(x);
  1721.     if (f)
  1722.     {
  1723.       avma=l;
  1724.       p1=cgetg(3,11);p1[1]=0x7ffe +lx;
  1725.       *px=p1;
  1726.     }
  1727.     else
  1728.     {
  1729.       tetpil=avma;
  1730.       p1=cgetg(lx-i+2,11);
  1731.       p1[1]=0x1007ffe +e+i;
  1732.       for (j=i;j<lx;j++)
  1733.         p1[j-i+2]=lcopy(x[j]);
  1734.       *px=gerepile(l,tetpil,p1);
  1735.     }
  1736.     setvarn(*px,v);
  1737.   }
  1738.   else
  1739.   {
  1740.     if(f) setsigne(x,0);else setsigne(x,1);
  1741.   }
  1742. }
  1743.  
  1744. void normalizepol(px)
  1745.      
  1746.      GEN  *px;
  1747.      
  1748. {
  1749.   long    i,lx,f;
  1750.   GEN     x;
  1751.   
  1752.   if(typ(x= *px)!=10) err(gnormaler);
  1753.   lx=lgef(x);
  1754.   if(lx>2)
  1755.     {
  1756.       f=1;i=lx-1;
  1757.       while (f&&(i>1))
  1758.     {f=isexactzero(x[i]);i--;}
  1759.       if(f) {setlgef(*px,2);setsigne(*px,0);}
  1760.       else
  1761.     {
  1762.       i++;f=gcmp0(x[i]);setlgef(*px,i+1);
  1763.       while(f&&(i>1))
  1764.         {f=gcmp0(x[i]);i--;}
  1765.       setsigne(*px,f?0:1);
  1766.     }
  1767.     }
  1768.   else setsigne(*px,0);
  1769. }
  1770.  
  1771. long gsigne(x)
  1772.      GEN x;
  1773. {
  1774.   switch(typ(x))
  1775.   {
  1776.     case 1:
  1777.     case 2: return signe(x);
  1778.     case 4:
  1779.     case 5: return (signe(x[2])>0) ? signe(x[1]) : -signe(x[1]);
  1780.     default: err(gsigner);
  1781.   }
  1782. }
  1783.  
  1784. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1785. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1786. /*                                                                 */
  1787. /*                            PUISSANCE                            */
  1788. /*                                                                 */
  1789. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1790. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1791.  
  1792. GEN     gpui(x,n,prec)
  1793.      
  1794.      GEN     x,n;
  1795.      long    prec;
  1796.      
  1797. {
  1798.   long    av1,av2,av3,lx,tx,i,j,tetpil;
  1799.   unsigned long m;
  1800.   GEN     p1,p2,y,z,alp;
  1801.   
  1802.   if (typ(n)==1)
  1803.   {
  1804.     if((lgef(n)<=3)&&cmpis(n,0x7fffffff)<0) y=gpuigs(x,itos(n),prec);
  1805.     else
  1806.     {
  1807.       z=x;av1=avma;
  1808.       if((typ(x)!=16)&&(typ(x)!=15)) y=gun;
  1809.       else
  1810.       {
  1811.     if(typ(x)==16)
  1812.       {
  1813.         p1=mulii(x[1],x[3]);p2=shifti(mulii(x[2],x[2]),-2);
  1814.         y=cgetg(4,16);y[1]=un;y[2]=mpodd(x[2]) ? un : zero;
  1815.         y[3]=lsubii(p1,p2);
  1816.       }
  1817.     else
  1818.       {
  1819.         y=cgetg(5,15);y[1]=un;
  1820.         p1=subii(mulii(x[2],x[2]),shifti(mulii(x[1],x[3]),2));
  1821.         p2=racine(p1);if(mpodd(subii(p2,x[2]))) p2=addsi(-1,p2);
  1822.         y[2]=(long)p2;y[3]=lshifti(subii(mulii(p2,p2),p1),-2);
  1823.         p1=cgetr(prec);y[4]=(long)p1;p1[1]=0x800000-((prec-2)<<5);
  1824.         p1[2]=0;tetpil=avma;y=gerepile(av1,tetpil,gcopy(y));
  1825.       }
  1826.       }
  1827.       for (i=lgef(n)-1;i>2;i--)
  1828.       {
  1829.         for (m=n[i],j=0;j<32;j++,m>>=1)
  1830.         {
  1831.           if (m&1) y=gmul(y,z);
  1832.           z=gsqr(z);
  1833.         }
  1834.       }
  1835.       for (m=n[2];m>1;m>>=1)
  1836.       {
  1837.         if (m&1) y=gmul(y,z);
  1838.         z=gsqr(z);
  1839.       }
  1840.       if (signe(n)>0)
  1841.       {
  1842.         tetpil=avma;y=gmul(y,z);
  1843.       }
  1844.       else
  1845.       {
  1846.         y=gmul(y,z);
  1847.         tetpil=avma;y=ginv(y);
  1848.       }
  1849.       y=gerepile(av1,tetpil,y);
  1850.     }
  1851.   }
  1852.   else
  1853.   {
  1854.     if((tx=typ(x))>=17)
  1855.       {
  1856.     y=cgetg(lx=lg(x),tx);for(i=1;i<lx;i++) y[i]=lpui(x[i],n,prec);
  1857.       }
  1858.     else
  1859.       {
  1860.     if(tx!=11)
  1861.       {
  1862.         av1=avma;
  1863.         if(gcmp0(x))
  1864.           {
  1865.         if(gcmpgs(p1=greal(n),0)<=0) err(gpuier3);
  1866.         av2=precision(x);
  1867.         if(av2)
  1868.           {
  1869.             p1=ground(gmulsg(gexpo(x),p1));
  1870.             if(lgef(p1)>3) err(gpuier4);
  1871.             av2=itos(p1);if(abs(av2)>=0x800000) err(gpuier4);
  1872.             avma=av1;y=cgetr(3);y[2]=0;y[1]=0x800000+av2;
  1873.           }
  1874.         else {avma=av1;y=gcopy(x);}
  1875.           }
  1876.         else
  1877.           {
  1878.         y=gmul(n,glog(x,prec));av2=avma;
  1879.         y=gerepile(av1,av2,gexp(y,prec));
  1880.           }
  1881.       }
  1882.     else
  1883.       {
  1884.         if(valp(x)) err(gpuier2);
  1885.         if(gvar9(n)>varn(x))
  1886.           {
  1887.         if(gcmp1(x[2]))
  1888.           {
  1889.             av1=avma;alp=gaddgs(n,1);
  1890.             av2=avma;y=cgetg(lx=lg(x),11);y[1]=0x1008000;
  1891.             y[2]=un;setvarn(y,varn(x));
  1892.             for(i=3;i<lx;i++)
  1893.               {
  1894.             av3=avma;p1=gzero;
  1895.             for(j=1;j<i-1;j++)
  1896.               {
  1897.                 p2=gsubgs(gmulgs(alp,j),i-2);
  1898.                 p1=gadd(p1,gmul(gmul(p2,x[j+2]),y[i-j]));
  1899.               }
  1900.             tetpil=avma;
  1901.             y[i]=lpile(av3,tetpil,gdivgs(p1,i-2));
  1902.               }
  1903.             y=gerepile(av1,av2,y);
  1904.           }
  1905.         else
  1906.           {
  1907.             av1=avma;p1=gdiv(x,x[2]);
  1908.             p1=gpui(p1,n,prec);p2=gpui(x[2],n,prec);
  1909.             tetpil=avma;y=gerepile(av1,tetpil,gmul(p1,p2));
  1910.           }
  1911.           }
  1912.         else
  1913.           {
  1914.         y=gmul(n,glog(x,prec));av2=avma;
  1915.         y=gerepile(av1,av2,gexp(y,prec));
  1916.           }
  1917.       }
  1918.       }
  1919.   }
  1920.   return y;
  1921. }
  1922.  
  1923. GEN     gpuigs(x,n,prec)
  1924.      
  1925.      GEN     x;
  1926.      long    n,prec;
  1927.      
  1928. {
  1929.   long    lx,av,m,dd,tetpil,i,j;
  1930.   GEN     y,z,p1,p2;
  1931.   
  1932.   
  1933.   if (!n)
  1934.   {
  1935.     lx=lg(x);
  1936.     switch(typ(x))
  1937.     {
  1938.       case 1 :
  1939.       case 2 :
  1940.       case 4 :
  1941.       case 5 : y=gun;break;
  1942.       case 6 : y=cgetg(3,6);y[1]=un;
  1943.         y[2]=zero;
  1944.         break;
  1945.       case 3 : y=cgetg(3,3);y[1]=copyifstack(x[1]);
  1946.         y[2]=un;
  1947.         break;
  1948.       case 8 : y=cgetg(4,8);y[1]=copyifstack(x[1]);
  1949.         y[2]=un;y[3]=zero;
  1950.         break;
  1951.       case 10: 
  1952.       case 11: 
  1953.       case 13: 
  1954.       case 14: y=polun[varn(x)]; break;
  1955.       case 9: y=cgetg(3,9);y[1]=copyifstack(x[1]);
  1956.         y[2]=lpuigs(x[2],0);
  1957.         break;
  1958.       case 15 : av=avma;y=cgetg(5,15);y[1]=un;
  1959.         p1=subii(mulii(x[2],x[2]),shifti(mulii(x[1],x[3]),2));
  1960.     p2=racine(p1);if(mpodd(subii(p2,x[2]))) p2=addsi(-1,p2);
  1961.     y[2]=(long)p2;y[3]=lshifti(subii(mulii(p2,p2),p1),-2);
  1962.         p1=cgetr(prec);y[4]=(long)p1;p1[1]=0x800000-((prec-2)<<5);
  1963.     p1[2]=0;tetpil=avma;y=gerepile(av,tetpil,gcopy(y));break;
  1964.       case 16 : y=cgetg(4,16);y[1]=un;y[2]=mpodd(x[2]) ? un : zero;
  1965.         av=avma;p1=mulii(x[1],x[3]);p2=shifti(mulii(x[2],x[2]),-2);tetpil=avma;
  1966.         y[3]=lpile(av,tetpil,subii(p1,p2));break;
  1967.       case 17:
  1968.       case 18: err(gpuier1);
  1969.       case 19: if (lx!=(lg(x[1]))) err(gpuier1);
  1970.         else
  1971.         {
  1972.           y=cgetg(lx,19);
  1973.           for (j=1;j<lx;j++)
  1974.           {
  1975.             y[j]=lgetg(lx,18);
  1976.             for (i=1;i<lx;i++)
  1977.               coeff(y,i,j)=(i!=j) ? zero :
  1978.                 lpuigs(coeff(x,i,i),0);
  1979.           }
  1980.         }
  1981.     }
  1982.   }
  1983.   else if (n==1) y=gcopy(x);
  1984.   else if (n== -1) y=ginv(x);
  1985.   else
  1986.   {
  1987.     m=abs(n);
  1988.     if(ismonome(x))
  1989.     {
  1990.       j=lgef(x);
  1991.       dd=(j-3)*m+3;
  1992.       av=avma;p1=cgetg(dd,10);p1[1]=0x1000000+dd;setvarn(p1,varn(x));
  1993.       p1[dd-1]=lpuigs(x[j-1],m);
  1994.       for(i=2;i<dd-1;i++) p1[i]=zero;
  1995.       if(n<0)
  1996.       {
  1997.         tetpil=avma;y=cgetg(3,13);y[1]=(long)denom(p1[dd-1]);
  1998.         y[2]=lmul(p1,y[1]);y=gerepile(av,tetpil,y);
  1999.       }
  2000.       else y=p1;
  2001.     }
  2002.     else
  2003.     {
  2004.       z=x;av=avma;
  2005.       if((typ(x)!=16)&&(typ(x)!=15)) y=gun;
  2006.       else
  2007.       {
  2008.     if(typ(x)==16)
  2009.       {
  2010.         p1=mulii(x[1],x[3]);p2=shifti(mulii(x[2],x[2]),-2);
  2011.         y=cgetg(4,16);y[1]=un;if(mpodd(x[2])) y[2]=un;else y[2]=zero;
  2012.         y[3]=lsubii(p1,p2);
  2013.       }
  2014.     else
  2015.       {
  2016.         y=cgetg(5,15);y[1]=un;
  2017.         p1=subii(mulii(x[2],x[2]),shifti(mulii(x[1],x[3]),2));
  2018.         p2=racine(p1);if(mpodd(subii(p2,x[2]))) p2=addsi(-1,p2);
  2019.         y[2]=(long)p2;y[3]=lshifti(subii(mulii(p2,p2),p1),-2);
  2020.         p1=cgetr(prec);y[4]=(long)p1;p1[1]=0x800000-((prec-2)<<5);
  2021.         p1[2]=0;tetpil=avma;y=gerepile(av,tetpil,gcopy(y));
  2022.       }
  2023.       }
  2024.       for(;m>1;m>>=1)
  2025.       {
  2026.         if (m&1) y=gmul(y,z);
  2027.         z=gsqr(z);
  2028.       }
  2029.       if (n>0)
  2030.       {
  2031.         tetpil=avma;y=gmul(y,z);
  2032.       }
  2033.       else
  2034.       {
  2035.         y=gmul(y,z);
  2036.         tetpil=avma;y=ginv(y);
  2037.       }
  2038.       y=gerepile(av,tetpil,y);
  2039.     }
  2040.   }
  2041.   return y;
  2042. }
  2043.  
  2044. long iscomplex(x)
  2045.   GEN x;
  2046.   
  2047. {
  2048.   long tx=typ(x);
  2049.   
  2050.   switch(tx)
  2051.   {
  2052.     case 1:
  2053.     case 2:
  2054.     case 4:
  2055.     case 5: return 0;
  2056.     case 6: return !gcmp0(x[2]);
  2057.     case 8: return signe(((GEN)x[1])[2])>0;
  2058.     default: err(iscomplexer1);
  2059.   }
  2060. }
  2061.  
  2062. long ismonome(x)
  2063.      GEN x;
  2064.  
  2065. {
  2066.   long i,l;
  2067.  
  2068.   if((typ(x)!=10)||(!signe(x))) return 0;
  2069.   l=lgef(x)-1;i=2;
  2070.   while((i<=l)&&isexactzero(x[i])) i++;
  2071.   return i==l;
  2072. }
  2073.  
  2074. GEN     gsqr(x)
  2075.      
  2076.      GEN     x;
  2077.      
  2078. {
  2079.   long  tx=typ(x),lx,av,i,j;
  2080.   GEN   y,p1;
  2081.   
  2082.   if(tx==7)
  2083.   {
  2084.     y=cgetg(5,7);
  2085.     y[2]=x[2];
  2086.     if(!cmpis(x[2] ,2))
  2087.     {
  2088.       i=precp(x)+1;av=avma;
  2089.       p1=addii(x[3],shifti(x[4],1));
  2090.       if(!gcmp0(p1))
  2091.       {
  2092.         j=vali(p1);if(j<i) i=j;
  2093.       }
  2094.       avma=av;y[3]=lshifti(x[3],i);
  2095.       setprecp(y,precp(x)+i);
  2096.     }
  2097.     else
  2098.     {
  2099.       y[3]=lcopy(x[3]);
  2100.       setprecp(y,precp(x));
  2101.     }
  2102.     y[4]=lgeti(lg(y[3]));
  2103.     setvalp(y,2*valp(x));av=avma;
  2104.     modiiz(mulii(x[4],x[4]),y[3],y[4]);
  2105.     avma=av;
  2106.   }
  2107.   else if(tx==15) y=sqcompreal(x);
  2108.   else if(tx==16) y=sqcomp(x);
  2109.   else
  2110.   {
  2111.     if((tx<17)||((tx==19)&&(lg(x)==lg(x[1])))) y=gmul(x,x);
  2112.     else
  2113.     {
  2114.       y=cgetg(lx=lg(x),tx);
  2115.       for(i=1;i<lx;i++)
  2116.         y[i]=lsqr(x[i]);
  2117.     }
  2118.   }
  2119.   return y;
  2120. }
  2121.