home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 2 / RISC_DISC_2.iso / pd_share / program / code / pari / c / alglin2 < prev    next >
Encoding:
Text File  |  1991-12-08  |  33.4 KB  |  1,362 lines
  1. /********************************************************************/
  2. /********************************************************************/
  3. /**                                                                **/
  4. /**                  ++++++++++++++++++++++++++++++                **/
  5. /**                  +                            +                **/
  6. /**                  +     ALGEBRE LINEAIRE       +                **/
  7. /**                  +                            +                **/
  8. /**                  ++++++++++++++++++++++++++++++                **/
  9. /**                                                                **/
  10. /**                        (deuxieme partie)                       **/
  11. /**                                                                **/
  12. /**                       copyright Babe Cool                      **/
  13. /**                                                                **/
  14. /**                                                                **/
  15. /********************************************************************/
  16. /********************************************************************/
  17.  
  18. # include "genpari.h"
  19.  
  20. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  21. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  22. /*                                                                 */
  23. /*                      POLYNOME CARACTERISTIQUE                   */
  24. /*                                                                 */
  25. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  26. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  27.  
  28. GEN     caract(x,v)
  29.      
  30.      GEN     x;
  31.      int     v;
  32.      
  33. {
  34.   long    n,k,l,tetpil,tx=typ(x);
  35.   GEN     p1,p2,p3,p4,p5,p6;
  36.   
  37.   switch(tx)
  38.     {
  39.     case 1: case 2: case 3: case 4: case 5: case 7: p1=cgetg(4,10);
  40.       p1[1]=0x1000004+(v<<16);p1[2]=lneg(x);p1[3]=un;return p1;
  41.     case 6: case 8: p1=cgetg(5,10);p1[1]=0x1000005+(v<<16);p1[2]=lnorm(x);
  42.       l=avma;p2=trace(x[1]);tetpil=avma;p1[3]=lpile(l,tetpil,gneg(p2));
  43.       p1[4]=un;return p1;
  44.     case 9: l=avma;p1=gsub(polx[255],x[2]);tetpil=avma;
  45.       p1=subres(x[1],p1);
  46.       if(varn(p1)==255) {setvarn(p1,v);return gerepile(l,tetpil,p1);}
  47.       else {tetpil=avma;return gerepile(l,tetpil,gsubst(p1,255,polx[v]));}
  48.     case 19: n=lg(x)-1;if(!n) return polun[v];
  49.       if((n+1)!=lg(x[1])) err(mattype1);
  50.       l=avma;p1=gzero;p2=gun;
  51.       if(n%2) p2=gneg(p2);p5=cgetg(4,10);
  52.       p5[1]=0x01000004;p5[3]=un;setvarn(p5,v);
  53.       p6=cgeti(3);p5[2]=zero;
  54.       p6[1]=0xff000003;p4=cgetg(3,14);p4[2]=(long)p5;
  55.       for(k=0;k<=n;k++)
  56.     {
  57.       p3=det(gsub(gscalsmat(k,n),x));p4[1]=lmul(p3,p2);p6[2]=k;
  58.       p1=gadd(p4,p1);p5[2]=(long)p6;
  59.       if(k!=n) p2=gdivgs(gmulsg(k-n,p2),k+1);
  60.     }
  61.       p2=mpfact(n);tetpil=avma;return gerepile(l,tetpil,gdiv(p1[1],p2));
  62.     default: err(mattype1);
  63.     }
  64. }
  65.  
  66.  
  67. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  68. /*                      Methode des traces :
  69.                         ce programme retourne le polynome caracteristique,
  70.                         et si un pointeur non nul est fourni,celui pointe
  71.                         sur la matrice adjointe a la sortie.       */
  72. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  73.  
  74. GEN     caradj(x,v,py)
  75.      
  76.      GEN     x,*py;
  77.      long    v;
  78.      
  79. {
  80.   long  i,j,k,l,t,av1,av2,av3,av4,decal;
  81.   GEN   p,y,z;
  82.   
  83.   if(typ(x)!=19) err(mattype1);
  84.   l=lg(x);
  85.   p=cgetg(l+2,10);setvarn(p,v);
  86.   p[l+1]=un;
  87.   av1=avma;t=ltrace(x);av2=avma;
  88.   t=lpile(av1,av2,gneg(t));
  89.   p[l]=t;
  90.   av1=avma;
  91.   y=cgetg(l,19);
  92.   for (i=1;i<l;i++) y[i]=lgetg(l,18);
  93.   for (i=1;i<l;i++)
  94.     for (j=1;j<l;j++)
  95.       {
  96.         if (i==j) coeff(y,i,j)=ladd(coeff(x,i,j),t);
  97.         else coeff(y,i,j)=coeff(x,i,j);
  98.       }
  99.   
  100.   for (k=2;k<l-1;k++)
  101.     {
  102.       z=gmul(x,y);
  103.       t=ltrace(z);av2=avma;
  104.       t=ldivgs(t,-k);av3=avma;
  105.       y=cgetg(l,19);
  106.       for (i=1;i<l;i++) y[i]=lgetg(l,18);
  107.       for (i=1;i<l;i++) for (j=1;j<l;j++)
  108.         if (i==j) coeff(y,i,j)=ladd(coeff(z,i,j),t);
  109.         else coeff(y,i,j)=lcopy(coeff(z,i,j));
  110.       av4=avma;decal=lpile(av1,av2,0);
  111.       p[l-k+1]=adecaler(t,av2,av4)?t+decal:t;
  112.       if(adecaler(y,av2,av4)) y+=(decal>>2);
  113.       av1=av3+decal;
  114.     }
  115.   t=zero;
  116.   for (i=1;i<l;i++)
  117.     t=ladd(t,gmul(coeff(x,1,i),coeff(y,i,1)));
  118.   av2=avma;t=lneg(t);
  119.   if ((long) py)
  120.     {
  121.       z=cgetg(l,19);
  122.       for (i=1;i<l;i++) z[i]=lgetg(l,18);
  123.       for (i=1;i<l;i++) for (j=1;j<l;j++)
  124.         coeff(z,i,j)=lcopy(coeff(y,i,j));
  125.       av4=avma;decal=lpile(av1,av2,0);
  126.       p[2]=adecaler(t,av2,av4)?t+decal:t;
  127.       *py=adecaler(z,av2,av4)?z+(decal>>2):z;
  128.     }
  129.   else p[2]=lpile(av1,av2,t);
  130.   p[1]=0x01000000+l+2;
  131.   return p;
  132. }
  133.  
  134.  
  135. GEN     adj(x)
  136.      GEN     x;
  137.      
  138. {
  139.   GEN     y;
  140.   
  141.   caradj(x,255,&y);
  142.   return y;
  143. }
  144.  
  145. GEN  caradj0(x,v)
  146.      GEN x;
  147.      long v;
  148. {
  149.   long tx=typ(x),l,tetpil;
  150.   GEN p1,p2;
  151.  
  152.   switch(tx)
  153.     {
  154.     case 1: case 2: case 3: case 4: case 5: case 7: p1=cgetg(4,10);
  155.       p1[1]=0x1000004+(v<<16);p1[2]=lneg(x);p1[3]=un;return p1;
  156.     case 6: case 8: p1=cgetg(5,10);p1[1]=0x1000005+(v<<16);p1[2]=lnorm(x);
  157.       l=avma;p2=trace(x);tetpil=avma;p1[3]=lpile(l,tetpil,gneg(p2));
  158.       p1[4]=un;return p1;
  159.     case 9: l=avma;p1=gsub(polx[255],x[2]);tetpil=avma;
  160.       p1=subres(x[1],p1);
  161.       if(varn(p1)==255) {setvarn(p1,v);return gerepile(l,tetpil,p1);}
  162.       else {tetpil=avma;return gerepile(l,tetpil,gsubst(p1,255,polx[v]));}
  163.     default: return caradj(x,v,0);
  164.     }
  165. }
  166.  
  167. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  168. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  169. /*                                                                 */
  170. /*              INVERSION D'UNE MATRICE                            */
  171. /*                                                                 */
  172. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  173. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  174.  
  175.  
  176. GEN     invmulmat(a,b)     /* calcule a^(-1)*b, b etant une matrice.
  177.                   ( Il faut : nblig(a)=nbcol(a)=nblig(b) ) */
  178.      GEN     a,b;
  179.      
  180. {
  181.   long  p,u,m,nbli,nbco,i,j,k,av,av1,av2,av3,av4;
  182.   GEN   aa,x;
  183.   
  184.   
  185.   nbco=lg(b)-1;if(!nbco) return cgetg(1,19);
  186.   nbli=lg(a[1])-1;
  187.   if (lg(a)-1 != nbli) err(invmuler1);
  188.   if (nbli!=lg(b[1])-1) err(invmuler1);
  189.   av=avma;
  190.   x=cgetg(nbco+1,19);
  191.   for (j=1;j<=nbco;j++)
  192.     {
  193.       x[j]=lgetg(nbli+1,18);
  194.       for (i=1;i<=nbli;i++)
  195.     coeff(x,i,j)=coeff(b,i,j);
  196.     }
  197.   av1=avma;
  198.   aa=cgetg(nbli+1,19);
  199.   for (j=1;j<=nbli;j++)
  200.     {
  201.       aa[j]=lgetg(nbli+1,18);
  202.       for (i=1;i<=nbli;i++) coeff(aa,i,j)=coeff(a,i,j);
  203.     }
  204.   for (i=1;i<nbli;i++)
  205.     {
  206.       p=coeff(aa,i,i);k=i;
  207.       if (gcmp0(p))
  208.     {
  209.       for (k=i+1;(k<=nbli)&&gcmp0(coeff(aa,k,i));k++);
  210.       if (k>nbli) err(matinv2);
  211.       else
  212.         {
  213.           for (j=i;j<=nbli;j++)
  214.         {
  215.           u=coeff(aa,i,j);coeff(aa,i,j)=coeff(aa,k,j);
  216.           coeff(aa,k,j)=u;
  217.         }
  218.           for (j=1;j<=nbco;j++)
  219.         {
  220.           u=coeff(x,i,j);coeff(x,i,j)=coeff(x,k,j);
  221.           coeff(x,k,j)=u;
  222.         }
  223.           p=coeff(aa,i,i);
  224.         }
  225.     }
  226.       for (k=i+1;k<=nbli;k++)
  227.     {
  228.       m=coeff(aa,k,i);
  229.       if (!gcmp0(m))
  230.         {
  231.           m=ldiv(m,p);
  232.           for (j=i+1;j<=nbli;j++)
  233.         coeff(aa,k,j)=lsub(coeff(aa,k,j),gmul(m,coeff(aa,i,j)));
  234.           for (j=1;j<=nbco;j++)
  235.         coeff(x,k,j)=lsub(coeff(x,k,j),gmul(m,coeff(x,i,j)));
  236.         }
  237.     }
  238.     }
  239.   av2=avma;
  240.   p=coeff(aa,nbli,nbli);
  241.   if (gcmp0(p)) err(matinv2);
  242.   else
  243.     {
  244.       for (j=1;j<=nbco;j++)
  245.     {
  246.       coeff(x,nbli,j)=ldiv(coeff(x,nbli,j),p);
  247.       
  248.       for (i=nbli-1;i>0;i--)
  249.         {
  250.           av3=avma;
  251.           m=coeff(x,i,j);
  252.           for (k=i+1;k<=nbli;k++)
  253.         m= lsub(m,gmul(coeff(aa,i,k),coeff(x,k,j)));
  254.           av4=avma;
  255.           coeff(x,i,j)=lpile(av3,av4,gdiv(m,coeff(aa,i,i)));
  256.         }
  257.     }
  258.       av=lpile(av1,av2,0);
  259.       for(i=1;i<=nbli;i++)
  260.     for(j=1;j<=nbco;j++)
  261.       if (coeff(x,i,j)<=av1) coeff(x,i,j)+=av;
  262.     }
  263.   return x;
  264. }
  265.  
  266. GEN invmat(a)
  267.      GEN a;
  268.      
  269. {
  270.   long av=avma,tetpil;
  271.   GEN b;
  272.   
  273.   b=gscalmat(un,lg(a)-1);tetpil=avma;
  274.   return gerepile(av,tetpil,invmulmat(a,b));
  275. }
  276.  
  277.  
  278. GEN     invmulmatreel(a,b)     /* calcule a^(-1)*b, b etant une matrice.
  279.                   ( Il faut : nblig(a)=nbcol(a)=nblig(b) ) */
  280.      GEN     a,b;
  281.      
  282. {
  283.   long  p,u,m,nbli,nbco,i,j,k,av,av1,av2,av3,av4;
  284.   GEN   aa,x,p1;
  285.   
  286.   
  287.   nbco=lg(b)-1;nbli=lg(a[1])-1;
  288.   if (lg(a)-1 != nbli) err(invmuler1);
  289.   /* la verif nblig(b)=nblig(a) n'est pas faite ... */
  290.   av=avma;
  291.   x=cgetg(nbco+1,19);
  292.   for (j=1;j<=nbco;j++)
  293.     {
  294.       x[j]=lgetg(nbli+1,18);
  295.       for (i=1;i<=nbli;i++)
  296.     coeff(x,i,j)=coeff(b,i,j);
  297.     }
  298.   av1=avma;
  299.   aa=cgetg(nbli+1,19);
  300.   for (j=1;j<=nbli;j++)
  301.     {
  302.       aa[j]=lgetg(nbli+1,18);
  303.       for (i=1;i<=nbli;i++) coeff(aa,i,j)=coeff(a,i,j);
  304.     }
  305.   for (i=1;i<nbli;i++)
  306.     {
  307.       p=labs(coeff(aa,i,i));k=i;
  308.       for(j=i+1;j<=nbli;j++)
  309.     if(gcmp(p1=gabs(coeff(aa,j,i),3),p)>0) {p=(long)p1;k=j;}
  310.       if (gcmp0(p)) err(matinv2);
  311.       else
  312.     {
  313.       if(k>i)
  314.         {
  315.           for (j=i;j<=nbli;j++)
  316.         {
  317.           u=coeff(aa,i,j);coeff(aa,i,j)=coeff(aa,k,j);
  318.           coeff(aa,k,j)=u;
  319.         }
  320.           for (j=1;j<=nbco;j++)
  321.         {
  322.           u=coeff(x,i,j);coeff(x,i,j)=coeff(x,k,j);
  323.           coeff(x,k,j)=u;
  324.         }
  325.         }
  326.       p=coeff(aa,i,i);
  327.     }
  328.       for (k=i+1;k<=nbli;k++)
  329.     {
  330.       m=coeff(aa,k,i);
  331.       if (!gcmp0(m))
  332.         {
  333.           m=ldiv(m,p);
  334.           for (j=i+1;j<=nbli;j++)
  335.         coeff(aa,k,j)=lsub(coeff(aa,k,j),gmul(m,coeff(aa,i,j)));
  336.           for (j=1;j<=nbco;j++)
  337.         coeff(x,k,j)=lsub(coeff(x,k,j),gmul(m,coeff(x,i,j)));
  338.         }
  339.     }
  340.     }
  341.   av2=avma;
  342.   p=coeff(aa,nbli,nbli);
  343.   if (gcmp0(p)) err(matinv2);
  344.   else
  345.     {
  346.       for (j=1;j<=nbco;j++)
  347.     {
  348.       coeff(x,nbli,j)=ldiv(coeff(x,nbli,j),p);
  349.       
  350.       for (i=nbli-1;i>0;i--)
  351.         {
  352.           av3=avma;
  353.           m=coeff(x,i,j);
  354.           for (k=i+1;k<=nbli;k++)
  355.         m= lsub(m,gmul(coeff(aa,i,k),coeff(x,k,j)));
  356.           av4=avma;
  357.           coeff(x,i,j)=lpile(av3,av4,gdiv(m,coeff(aa,i,i)));
  358.         }
  359.     }
  360.       av=lpile(av1,av2,0);
  361.       for(i=1;i<=nbli;i++)
  362.     for(j=1;j<=nbco;j++)
  363.       if (coeff(x,i,j)<=av1) coeff(x,i,j)+=av;
  364.     }
  365.   return x;
  366. }
  367.  
  368. GEN invmatreel(a)
  369.      GEN a;
  370.      
  371. {
  372.   long av=avma,tetpil;
  373.   GEN b;
  374.   
  375.   b=gscalmat(un,lg(a)-1);tetpil=avma;
  376.   return gerepile(av,tetpil,invmulmatreel(a,b));
  377. }
  378.  
  379. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  380. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  381. /*                                                                 */
  382. /*              FORME DE HESSENBERG D'UNE MATRICE                  */
  383. /*                                                                 */
  384. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  385. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  386.  
  387. GEN hess(x)
  388.      GEN x;
  389. {
  390.   long tx=typ(x),lx=lg(x),av=avma,tetpil,m,i,j;
  391.   GEN p1,p2,p3,y;
  392.   
  393.   if((tx!=19)||(lg(x[1])!=lx)) err(mattype1);
  394.   y=gcopy(x);
  395.   for(m=2;m<lx-1;m++)
  396.     {
  397.       p2=gzero;
  398.       for(i=m+1;(i<lx)&&(gcmp0(p2=(GEN)coeff(y,i,m-1)));i++);
  399.       if(!gcmp0(p2))
  400.     {
  401.       for(j=m-1;j<lx;j++)
  402.         {
  403.           p1=(GEN)coeff(y,i,j);coeff(y,i,j)=coeff(y,m,j);coeff(y,m,j)=(long)p1;
  404.         }
  405.       p1=(GEN)y[i];y[i]=y[m];y[m]=(long)p1;
  406.       for(i=m+1;i<lx;i++)
  407.         {
  408.           if(!gcmp0(p3=(GEN)coeff(y,i,m-1)))
  409.         {
  410.           p3=gdiv(p3,p2);coeff(y,i,m-1)=zero;
  411.           for(j=m;j<lx;j++)
  412.             coeff(y,i,j)=lsub(coeff(y,i,j),gmul(p3,coeff(y,m,j)));
  413.           for(j=1;j<lx;j++)
  414.             coeff(y,j,m)=ladd(coeff(y,j,m),gmul(p3,coeff(y,j,i)));
  415.         }
  416.         }
  417.     }
  418.     }
  419.   tetpil=avma;return gerepile(av,tetpil,gcopy(y));
  420. }
  421.  
  422. GEN carhess(x,v)
  423.      GEN x;
  424.      long v;
  425. {
  426.   long av=avma,tetpil,tx=typ(x),lx=lg(x),r,i;
  427.   GEN *y,p1,p2,p3,p4;
  428.   
  429.   if((tx!=19)||(lg(x[1])!=lx)) err(mattype1);
  430.   y=(GEN *)newbloc(4*lx);
  431.   y[0]=polun[v];p1=hess(x);p2=polx[v];
  432.   for(r=1;r<lx;r++)
  433.     {
  434.       y[r]=gmul(y[r-1],gsub(p2,coeff(p1,r,r)));
  435.       p3=gun;p4=gzero;
  436.       for(i=1;i<r;i++)
  437.     {
  438.       p3=gmul(p3,coeff(p1,r-i+1,r-i));
  439.       p4=gadd(p4,gmul(gmul(p3,coeff(p1,r-i,r)),y[r-i-1]));
  440.     }
  441.       tetpil=avma;y[r]=gsub(y[r],p4);
  442.     }
  443.   p1=gerepile(av,tetpil,y[lx-1]);
  444.   killbloc(y);return p1;
  445. }
  446.  
  447.  
  448.  
  449. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  450. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  451. /*                                                                 */
  452. /*                          NORME                                  */
  453. /*                                                                 */
  454. /*            Cree un GEN pointant sur la norme de x               */
  455. /*                                                                 */
  456. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  457. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  458.  
  459. GEN     gnorm(x)
  460.      
  461.      GEN     x;
  462.      
  463. {
  464.   long    l,tx,lx,i,tetpil;
  465.   GEN     p1,p2,y;
  466.   
  467.   switch(tx=typ(x))
  468.     {
  469.     case 1 :
  470.     case 2 : y=mpmul(x,x);break;
  471.       
  472.     case 3 : err(normer1);
  473.       
  474.     case 4 :
  475.     case 5 : y=gmul(x,x);break;
  476.       
  477.     case 6 : l=avma;p1=gmul(x[1],x[1]);
  478.       p2=gmul(x[2],x[2]);
  479.       tetpil=avma;
  480.       y=gerepile(l,tetpil,gadd(p1,p2));break;
  481.       
  482.     case 8 : l=avma;p1=(GEN)x[1];
  483.       if (gcmp0(p1[3]))
  484.         {
  485.           p2=gmul(p1[2],gmul(x[3],x[3]));
  486.           p1=gmul(x[2],x[2]);
  487.           tetpil=avma;
  488.           y=gerepile (l,tetpil,gadd(p1,p2));
  489.         }
  490.       else
  491.         {
  492.           p2=gmul(p1[2],gmul(x[3],x[3]));
  493.           p1=gmul(x[2],gadd(x[2],x[3]));
  494.           tetpil=avma;
  495.           y=gerepile(l ,tetpil,gadd(p1,p2));
  496.         }
  497.       break;
  498.     case 10:
  499.     case 11:
  500.     case 13:
  501.     case 14: l=avma;p1=gmul(gconj(x),x);tetpil=avma;
  502.       y=gerepile(l,tetpil,greal(p1));break;
  503.     case 9 : y=subres(x[1],x[2]);break;
  504.     case 17: 
  505.     case 18:
  506.     case 19: lx=lg(x);y=cgetg(lx,tx);
  507.       for(i=1;i<lx;i++) y[i]=lnorm(x[i]);
  508.       break;
  509.     default: err(normer1);
  510.     }
  511.   return y;
  512. }
  513.  
  514. GEN     gnorml2(x)
  515.      GEN   x;
  516.      
  517. {
  518.   GEN  y,p1;
  519.   long i,tx=typ(x),lx=lg(x),av,tetpil;
  520.   
  521.   if(tx<17) return gnorm(x);
  522.   y=gzero;
  523.   if(lx>1)
  524.     {
  525.       av=avma;y=gnorml2(x[1]);
  526.       for(i=2;i<lx;i++)
  527.         {
  528.           p1=gnorml2(x[i]);tetpil=avma;
  529.           y=gadd(p1,y);
  530.         }
  531.       if(lx>2) y=gerepile(av,tetpil,y);
  532.     }
  533.   return y;
  534. }
  535.  
  536.  
  537. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  538. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  539. /*                                                                 */
  540. /*                            CONJUGAISON                          */
  541. /*                                                                 */
  542. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  543. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  544.  
  545.  
  546. GEN     gconj(x)
  547.      
  548.      GEN     x;
  549.      
  550. {
  551.   long    lx,i,tx=typ(x);
  552.   GEN     z,p1;
  553.   
  554.   switch(tx)
  555.     {
  556.     case 1 :
  557.     case 2 :
  558.     case 3 :
  559.     case 4 :
  560.     case 5 :
  561.     case 7 : z=gcopy(x);break;
  562.       
  563.     case 6 : z=cgetg(3,6);z[1]=lcopy(x[1]);
  564.       z[2]=lneg(x[2]);
  565.       break;
  566.       
  567.     case 8 : z=cgetg(4,8);z[1]=copyifstack(x[1]);
  568.       p1=(GEN)x[1];
  569.       z[3]=lneg(x[3]);
  570.       if(gcmp0(p1[3])) z[2]=lcopy(x[2]);
  571.       else z[2]=ladd(x[2],x[3]);
  572.       break;
  573.       
  574.     case 10: lx=lg(x);z=cgetg(lx,tx);
  575.       z[1]=x[1];
  576.       for(i=2;i<lgef(x);i++)
  577.         z[i]=lconj(x[i]);
  578.       break;
  579.       
  580.     case 11: lx=lg(x);z=cgetg(lx,tx);
  581.       z[1]=x[1];
  582.       if(!gcmp0(x))
  583.         {
  584.           for(i=2;i<lx;i++)
  585.             z[i]=lconj(x[i]);
  586.         }
  587.       break;
  588.       
  589.     case 13:
  590.     case 14: 
  591.     case 17:
  592.     case 18:
  593.     case 19: lx=lg(x);z=cgetg(lx,tx);
  594.       for(i=1;i<lx;i++)
  595.         z[i]=lconj(x[i]);
  596.       break;
  597.       
  598.     default: err(conjer1);
  599.     }
  600.   return z;
  601. }
  602.  
  603. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  604. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  605. /*                                                                 */
  606. /*                  PARTIES REELLE ET IMAGINAIRES                  */
  607. /*                                                                 */
  608. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  609. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  610.  
  611.  
  612. GEN     greal(x)
  613.      
  614.      GEN     x;
  615.      
  616. {
  617.   long    lx,i,j,av,tetpil,tx=typ(x);
  618.   GEN     p1,p2,z;
  619.   
  620.   switch(tx)
  621.     {
  622.     case 1 :
  623.     case 2 :
  624.     case 4 :
  625.     case 5 : z=gcopy(x);break;
  626.       
  627.     case 6 : z=gcopy(x[1]);break;
  628.       
  629.     case 8 : z=gcopy(x[2]);break;
  630.       
  631.     case 10: lx=lgef(x);av=avma;
  632.       for(i=lx-1;(i>=2)&&(gcmp0(greal(x[i])));i--);
  633.       avma=av;
  634.       if(i<2) {z=cgetg(2, tx);z[1]=2;}
  635.       else 
  636.     {
  637.       z=cgetg(i+1,tx);z[1]=0x01000001+i;
  638.       for(j=2;j<=i;j++) z[j]=lreal(x[j]);
  639.     }
  640.       setvarn(z,varn(x));
  641.       break;
  642.       
  643.     case 11: lx=lg(x);av=avma;
  644.       if(gcmp0(x)) {z=cgetg(2,tx);z[1]=x[1];}
  645.       else
  646.     {
  647.       for(i=2;(i<lx)&&(gcmp0(greal(x[i])));i++);
  648.       avma=av;
  649.       if(i==lx)
  650.         {
  651.           z=cgetg(2,tx); setvalp(z, lx - 2 + valp(x));
  652.           setsigne(z,0); setvarn(z, varn(x));
  653.         }
  654.       else
  655.         {
  656.           z=cgetg(lx-i+2,tx);
  657.           for(j = 2; j <= lx - i + 1; j++) z[j] = lreal(x[j + i - 2]);
  658.           z[1] = x[1]; setvalp(z, valp(x) + i - 2);
  659.         }
  660.     }
  661.       break;
  662.       
  663.     case 13:
  664.     case 14: av=avma;p1=gadd(gmul(greal(x[1]),greal(x[2])),gmul(gimag(x[1]),gimag(x[2])));
  665.       p2=gadd(gsqr(greal(x[2])),gsqr(gimag(x[2])));tetpil=avma;
  666.       z=gerepile(av,tetpil,gdiv(p1,p2));break;
  667.     case 17:
  668.     case 18:
  669.     case 19: lx=lg(x);z=cgetg(lx,tx);
  670.       for(i=1;i<lx;i++) z[i]=lreal(x[i]);
  671.       break;
  672.       
  673.     default: err(realer1);
  674.     }
  675.   return z;
  676. }
  677.  
  678. GEN     gimag(x)
  679.      
  680.      GEN     x;
  681.      
  682. {
  683.   long    lx,i,j,av,tetpil,tx=typ(x);
  684.   GEN     p1,p2,z;
  685.   
  686.   switch(tx)
  687.     {
  688.     case 1 :
  689.     case 2 :
  690.     case 4 :
  691.     case 5 : z=gzero;break;
  692.       
  693.     case 6 : z=gcopy(x[2]);
  694.       break;
  695.       
  696.     case 8 : z=gcopy(x[3]);
  697.       break;
  698.       
  699.     case 10: lx=lgef(x);av=avma;
  700.       for(i=lx-1;(i>=2)&&(gcmp0(gimag(x[i])));i--);
  701.       avma=av;
  702.       if(i<2) {z=cgetg(2, tx);z[1]=2;}
  703.       else 
  704.     {
  705.       z=cgetg(i+1,tx);z[1]=0x01000001+i;
  706.       for(j=2;j<=i;j++) z[j]=limag(x[j]);
  707.     }
  708.       setvarn(z,varn(x));
  709.       break;
  710.       
  711.     case 11: lx=lg(x);av=avma;
  712.       if(gcmp0(x)) {z=cgetg(2,tx);z[1]=x[1];}
  713.       else
  714.     {
  715.       for(i=2;(i<lx)&&(gcmp0(gimag(x[i])));i++);
  716.       avma=av;
  717.       if(i==lx)
  718.         {
  719.           z=cgetg(2,tx); setvalp(z, lx - 2 + valp(x));
  720.           setsigne(z,0); setvarn(z, varn(x));
  721.         }
  722.       else
  723.         {
  724.           z=cgetg(lx-i+2,tx);
  725.           for(j = 2; j <= lx - i + 1; j++) z[j] = limag(x[j + i - 2]);
  726.           z[1] = x[1]; setvalp(z, valp(x) + i - 2);
  727.         }
  728.     }
  729.       break;
  730.       
  731.     case 13:
  732.     case 14: av=avma;p1=gsub(gmul(gimag(x[1]),greal(x[2])),gmul(greal(x[1]),gimag(x[2])));
  733.       p2=gadd(gsqr(greal(x[2])),gsqr(gimag(x[2])));tetpil=avma;
  734.       z=gerepile(av,tetpil,gdiv(p1,p2));break;
  735.     case 17:
  736.     case 18:
  737.     case 19: lx=lg(x);z=cgetg(lx,tx);
  738.       for(i=1;i<lx;i++)
  739.         z[i]=limag(x[i]);
  740.       break;
  741.       
  742.     default: err(imager1);
  743.     }
  744.   return z;
  745. }
  746.  
  747. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  748. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  749. /*                                                                 */
  750. /*                            TRACES                               */
  751. /*                                                                 */
  752. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  753. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  754.  
  755. GEN     assmat(x)
  756.      
  757.      GEN     x;
  758.      
  759. {
  760.   long    lx,i,j;
  761.   GEN     y,p1;
  762.   
  763.   if((typ(x)!=10) || gcmp0(x)) err(mattype2);
  764.   lx=lgef(x)-2;y=cgetg(lx,19);
  765.   for(i=1;i<lx-1;i++)
  766.     {
  767.       p1=cgetg(lx,18);y[i]=(long)p1;
  768.       for(j=1;j<lx;j++)
  769.         {p1[j]=(j==(i+1)) ? un : zero;}
  770.     }
  771.   p1=cgetg(lx,18);y[i]=(long)p1;
  772.   if(gcmp1(x[lx+1]))
  773.     {
  774.       for(j=1;j<lx;j++)
  775.         p1[j]=lneg(x[j+1]);
  776.     }
  777.   else
  778.     {
  779.       gnegz(x[lx+1],x[lx+1]);
  780.       for(j=1;j<lx;j++)
  781.         p1[j]=ldiv(x[j+1],x[lx+1]);
  782.       gnegz(x[lx+1],x[lx+1]);
  783.     }
  784.   return y;
  785. }
  786.  
  787. GEN     trace(x)
  788.      
  789.      GEN     x;
  790.      
  791. {
  792.   long    i,l,n,tx=typ(x),lx=lg(x),tetpil;
  793.   GEN     y,p1,p2;
  794.   
  795.   switch(tx)
  796.     {
  797.     case 1 :
  798.     case 2 :
  799.     case 4 :
  800.     case 5 : y=gmul2n(x,1);break;
  801.       
  802.     case 6 : y=gmul2n(x[1],1);break;
  803.       
  804.     case 8 : p1=(GEN)x[1];
  805.       if (!gcmp0(p1[3]))
  806.         {
  807.           l=avma;p2=gmul2n(x[2],1);
  808.           tetpil=avma;
  809.           y=gerepile(l,tetpil,gadd(x[3],p2));
  810.         }
  811.       else y=gmul2n(x[2],1);
  812.       break;
  813.       
  814.     case 10: lx=lg(x);y=cgetg(lx,tx);
  815.       y[1]=x[1];
  816.       for(i=2;i<lgef(x);i++)
  817.         y[i]=ltrace(x[i]);
  818.       break;
  819.       
  820.     case 11: lx=lg(x);y=cgetg(lx,tx);
  821.       y[1]=x[1];
  822.       if(!gcmp0(x))
  823.         {
  824.           for(i=2;i<lx;i++)
  825.             y[i]=ltrace(x[i]);
  826.         }
  827.       break;
  828.       
  829.     case 9 : l=avma;p1=polsym(x[1],n=(lgef(x[1])-4));
  830.       p2=gzero;for(i=0;i<=n;i++) p2=gadd(p2,gmul(truecoeff(x[2],i),p1[i+1]));
  831.       tetpil=avma;y=gerepile(l,tetpil,gcopy(p2));
  832.       break;
  833.       
  834.     case 13:
  835.     case 14: y=gadd(x,gconj(x));break;
  836.       
  837.     case 17:
  838.     case 18: lx=lg(x);y=cgetg(lx,tx);
  839.       for(i=1;i<lx;i++)
  840.         y[i]=ltrace(x[i]);
  841.       break;
  842.       
  843.     case 19: if (lx!=lg(x[1])) err(mattype3);
  844.       l=avma;p1=gcopy(coeff(x,1,1));
  845.       if(lx==2) return p1;
  846.       else
  847.         {
  848.           for(i=2;i<lx-1;i++)
  849.             p1=gadd(p1,coeff(x,i,i));
  850.           tetpil=avma;
  851.           y=gerepile(l,tetpil,gadd(p1,coeff(x,i,i)));
  852.         }
  853.       break;
  854.     default: err(mattype3);
  855.     }
  856.   return y;
  857. }
  858.  
  859. /*===================================*/
  860. /*     Reduction en carres           */
  861. /*===================================*/
  862.  
  863.  
  864. /*=======================================================
  865.   Reduction de Gauss ( Matrice definie >0 )
  866.   ========================================================*/
  867.  
  868.  
  869. GEN sqred1(a)
  870.      GEN a;
  871.      
  872. {
  873.   GEN b;
  874.   long av,av1,n,i,j,k,p,lim;
  875.   
  876.   if (typ(a)!=19) err(kerer1);
  877.   if (lg(a[1])!=(n=lg(a))) err(mattype1);
  878.   lim=(avma+bot)>>1;
  879.   n--;av=avma;
  880.   b=gcopy(a);
  881.   for(i=1;i<=n;i++)
  882.     for(j=1;j<i;j++) coeff(b,i,j) = zero;
  883.   for(k=1;k<=n;k++)
  884.     {
  885.       if(gsigne(p=coeff(b,k,k))<=0) err(sqreder1);
  886.       for(i=k+1;i<=n;i++)
  887.     for(j=i;j<=n;j++)
  888.       coeff(b,i,j)=lsub(coeff(b,i,j),gdiv(gmul(coeff(b,k,i),coeff(b,k,j)),p));
  889.       for(j=k+1;j<=n;j++)
  890.     coeff(b,k,j)=ldiv(coeff(b,k,j),p);
  891.       if(avma<lim) {av1=avma;b=gerepile(av,av1,gcopy(b));}
  892.     }
  893.   av1=avma;
  894.   return gerepile(av,av1,gcopy(b));
  895. }
  896.  
  897. GEN sqred3(a)
  898.      GEN a;
  899. {
  900.   long n,av=avma,tetpil,lim,i,j,k,l;
  901.   GEN p1,z;
  902.  
  903.   if (typ(a)!=19) err(kerer1);
  904.   if (lg(a[1])!=(n=lg(a))) err(mattype1);
  905.   lim=(avma+bot)>>1;
  906.   av=avma;z=cgetg(n,19);
  907.   for(j=1;j<n;j++) 
  908.     {
  909.       p1=cgetg(n,18);z[j]=(long)p1;
  910.       for(i=1;i<n;i++) p1[i]=zero;
  911.     }
  912.   for(i=1;i<n;i++)
  913.     {
  914.       for(k=1;k<i;k++)
  915.     {
  916.       p1=gzero;
  917.       for(l=1;l<k;l++) p1=gadd(p1,gmul(gmul(coeff(z,l,l),coeff(z,k,l)),coeff(z,i,l)));
  918.       coeff(z,i,k)=ldiv(gsub(coeff(a,i,k),p1),coeff(z,k,k));
  919.     }
  920.       p1=gzero;
  921.       for(l=1;l<i;l++) p1=gadd(p1,gmul(gmul(coeff(z,l,l),coeff(z,i,l)),coeff(z,i,l)));
  922.       coeff(z,i,k)=lsub(coeff(a,i,i),p1);
  923.       if(avma<lim) {tetpil=avma;z=gerepile(av,tetpil,gcopy(z));}
  924.     }
  925.   tetpil=avma;return gerepile(av,tetpil,gcopy(z));
  926. }
  927.  
  928. /*=======================================================
  929.   Reduction de Gauss (matrice symetrique quelconque)
  930.   Signature d'une matrice symetrique
  931.   ( seule la partie superieure est consideree )
  932.   ========================================================*/
  933.  
  934. GEN sqred2(a,flg)
  935.      GEN a; long flg;
  936.      
  937. {
  938.   GEN b,r;
  939.   long av,av1,av2,lim,n,i,j,k,l,p,sp,sn,t,u;
  940.   if (typ(a)!=19) err(kerer1);
  941.   if (lg(a[1])!=(n=lg(a))) err(mattype1);
  942.   av=avma;lim=(avma+bot)>>1;
  943.   r=cgeti(n);for(i=1;i<n;i++) r[i]=1;
  944.   av2=avma;b=gcopy(a);
  945.   t=(--n);sp=sn=0;
  946.   
  947.   while (t)
  948.     {
  949.       for(k=1;(k<=n)&&(gcmp0(coeff(b,k,k))||(!r[k]));k++);
  950.       if(k<=n)
  951.     {
  952.       p=coeff(b,k,k);
  953.       if(signe(p)>0) sp++;else sn++;
  954.       r[k]=0;t--;
  955.       for(j=1;j<=n;j++) 
  956.         coeff(b,k,j)=r[j] ? ldiv(coeff(b,k,j),p) : zero;
  957.       
  958.       for(i=1;i<=n;i++) if (r[i])
  959.         {
  960.           for(j=1;j<=n;j++)
  961.         coeff(b,i,j)=r[j] ? lsub(coeff(b,i,j),gmul(gmul(coeff(b,k,i),coeff(b,k,j)),p)) : zero;
  962.         }
  963.       coeff(b,k,k)=p;
  964.     }
  965.       else
  966.     {
  967.       for(k=1;k<=n;k++) if (r[k])
  968.         {
  969.           for(l=k+1;(l<=n)&&(gcmp0(coeff(b,k,l))||(!r[l]));l++);
  970.           if(l<=n)
  971.         {
  972.           p=coeff(b,k,l);r[k]=r[l]=0;sp++;sn++;t-=2;
  973.           for(i=1;i<=n;i++) if(r[i])
  974.             {
  975.               for(j=1;j<=n;j++)
  976.             coeff(b,i,j)=r[j]? lsub(coeff(b,i,j),gdiv(gadd(gmul(coeff(b,k,i),coeff(b,l,j)),gmul(coeff(b,k,j),coeff(b,l,i))),p)) : zero;
  977.             }
  978.           for(j=1;j<=n;j++) if (r[j])
  979.             {
  980.               u=coeff(b,k,j);
  981.               coeff(b,k,j)=ldiv(gadd(u,coeff(b,l,j)),p);
  982.               coeff(b,l,j)=ldiv(gsub(u,coeff(b,l,j)),p);
  983.             }
  984.           coeff(b,k,l)=un;coeff(b,l,k)=lneg(un);
  985.           coeff(b,k,k)=lmul2n(p,-1);coeff(b,l,l)=lneg(coeff(b,k,k));
  986.           break;
  987.         }
  988.           if(avma<lim) {av1=avma;b=gerepile(av2,av1,gcopy(b));}
  989.         }
  990.       if(k>n) break;
  991.     }
  992.     }
  993.   if (flg) {av1=avma;return gerepile(av,av1,gcopy(b));}
  994.   else
  995.     {
  996.       avma=av;
  997.       b=cgetg(3,17);b[1]=lstoi(sp);b[2]=lstoi(sn);return b;
  998.     }
  999. }
  1000.  
  1001. GEN sqred(a) { return sqred2(a,1); }
  1002. GEN signat(a) { return sqred2(a,0); }
  1003.  
  1004.  
  1005.  
  1006. /*===========================================================================
  1007.   Diagonalisation d'une matrice symetrique REELLE;
  1008.   matrice de passage orthogonale R
  1009.   (   Nouvelle version : 25/6/90  )
  1010.   (  Renvoie un vecteur a 2 comp :
  1011.   1-re comp = vect des valeurs propres
  1012.   2-me comp = matr des vecteurs propres   ).
  1013.   ============================================================================*/
  1014.  
  1015. GEN jacobi(a,prec) GEN a;long prec;
  1016.      
  1017. {
  1018.   long de,e,e1,e2,l,n,i,j,p,q,x,y,x1,y1,av1,av2,iter=0;
  1019.   GEN c,s,t,u,ja,lambda,r,unr;
  1020.   
  1021.   ja=cgetg(3,17);
  1022.   n=(l=lg(a))-1;
  1023.   e1=0x7fffff;
  1024.   lambda=cgetg(l,18);ja[1]=(long)lambda;
  1025.   for(j=1;j<=n;j++)
  1026.     {
  1027.       gaffect(coeff(a,j,j),x=lambda[j]=lgetr(prec));
  1028.       if(((e=expo(x))<e1)&&(gsigne(x))) e1=e;
  1029.     }
  1030.   r=cgetg(l,19);ja[2]=(long)r;
  1031.   for(j=1;j<=n;j++)
  1032.     {
  1033.       r[j]=lgetg(l,18);
  1034.       for(i=1;i<l;i++)
  1035.     affsr((i==j)? 1:0,coeff(r,i,j)=lgetr(prec));
  1036.     }
  1037.   av1=avma;
  1038.   e2=expo(coeff(a,1,2));p=1;q=2;
  1039.   
  1040.   c=cgetg(l,19);
  1041.   for(j=1;j<=n;j++)
  1042.     {
  1043.       c[j]=lgetg(j,18);
  1044.       for(i=1;i<j;i++)
  1045.     {
  1046.       gaffect(coeff(a,i,j),x=coeff(c,i,j)=lgetr(prec));
  1047.       if((e=expo(x))>e2) {e2=e;p=i;q=j;}
  1048.     }
  1049.     }
  1050.   
  1051.   a=c;
  1052.   affsr(1,unr=cgetr(prec));
  1053.   de=((prec-2)<<5);
  1054.   
  1055.   while(e1-e2<de)
  1056.     {
  1057.       iter++;
  1058.       /*calcul de la rotation associee dans le plan
  1059.     des p et q-iemes vecteurs de base   */
  1060.       av2=avma;
  1061.       x=ldivrr(subrr(lambda[q],lambda[p]),shiftr(coeff(a,p,q),1));
  1062.       y=lmpsqrt(addrr(unr,mulrr(x,x)));
  1063.       t=(gsigne(x)>0)? divrr(unr,addrr(x,y)) : divrr(unr,subrr(x,y));
  1064.       c=divrr(unr,mpsqrt(addrr(unr,mulrr(t,t))));
  1065.       s=mulrr(t,c);u=divrr(s,addrr(unr,c));
  1066.       
  1067.       /* Recalcul des transformees successives de la matrice a et de la matrice
  1068.      cumulee (r) des rotations :  */
  1069.       
  1070.       
  1071.       for(i=1;i<p;i++)
  1072.     {
  1073.       x=coeff(a,i,p); y=coeff(a,i,q);
  1074.       x1=lsubrr(x,mulrr(s,addrr(y,mulrr(u,x))));
  1075.       y1=laddrr(y,mulrr(s,subrr(x,mulrr(u,y))));
  1076.       affrr(x1,coeff(a,i,p));affrr(y1,coeff(a,i,q));
  1077.     }
  1078.       for(i=p+1;i<q;i++)
  1079.     {
  1080.       x=coeff(a,p,i); y=coeff(a,i,q);
  1081.       x1=lsubrr(x,mulrr(s,addrr(y,mulrr(u,x))));
  1082.       y1=laddrr(y,mulrr(s,subrr(x,mulrr(u,y))));
  1083.       affrr(x1,coeff(a,p,i));affrr(y1,coeff(a,i,q));
  1084.     }
  1085.       for(i=q+1;i<=n;i++)
  1086.     {
  1087.       x=coeff(a,p,i); y=coeff(a,q,i);
  1088.       x1=lsubrr(x,mulrr(s,addrr(y,mulrr(u,x))));
  1089.       y1=laddrr(y,mulrr(s,subrr(x,mulrr(u,y))));
  1090.       affrr(x1,coeff(a,p,i));affrr(y1,coeff(a,q,i));
  1091.     }
  1092.       x=lambda[p]; y=coeff(a,p,q); subrrz(x,mulrr(t,y),lambda[p]);
  1093.       x=y; y=lambda[q];  addrrz(y,mulrr(t,x),y);
  1094.       /*      if((e=expo(lambda[p]))<e1) e1=e;
  1095.           if((e=expo(lambda[q]))<e1) e1=e; */
  1096.       /*      affsr(0,x);      NON !  */
  1097.       setexpo(x,expo(x)-de-1);
  1098.       
  1099.       for(i=1;i<=n;i++)
  1100.     {
  1101.       x=coeff(r,i,p); y=coeff(r,i,q);
  1102.       x1=lsubrr(x,mulrr(s,addrr(y,mulrr(u,x))));
  1103.       y1=laddrr(y,mulrr(s,subrr(x,mulrr(u,y))));
  1104.       affrr(x1,coeff(r,i,p));affrr(y1,coeff(r,i,q));
  1105.     }
  1106.       
  1107.       e2=expo(coeff(a,1,2));p=1;q=2;
  1108.       for(j=1;j<=n;j++)
  1109.     {
  1110.       for(i=1;i<j;i++) if((e=expo(coeff(a,i,j)))>e2) {e2=e;p=i;q=j;}
  1111.       for(i=j+1;i<=n;i++) if((e=expo(coeff(a,j,i)))>e2) {e2=e;p=j;q=i;}
  1112.     }
  1113.       avma=av2;
  1114.     }     /* Fin de la boucle (while) de recalcul */
  1115.   avma=av1; return ja;
  1116. }
  1117.  
  1118. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1119. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1120. /*~                                    ~*/
  1121. /*~            MATRICE RATIONNELLE-->ENTIERE            ~*/
  1122. /*~                                    ~*/
  1123. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1124. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  1125.  
  1126. GEN matrixqz(x,pp)
  1127.      GEN x,pp;
  1128. {
  1129.   long av=avma,av1,tetpil,i,j,j1,m,n,t,fl,lim,nfact;
  1130.   GEN p,p1,p2,p3,p4,p5,unmodp,pk,pt;
  1131.  
  1132.   if(typ(x)!=19) err(matqzer1);
  1133.   lim=(avma+bot)>>1;
  1134.   n=lg(x)-1;if(!n) return gcopy(x);
  1135.   m=lg(x[1])-1;if(n>m) err(matqzer2);
  1136.   if(n==m)
  1137.     {
  1138.       p1=det(x);if(gcmp0(p1)) err(matqzer4);
  1139.       return idmat(n);
  1140.     }
  1141.   p1=cgetg(n+1,19);
  1142.   for(j=1;j<=n;j++)
  1143.     {
  1144.       p2=gun;p3=(GEN)x[j];
  1145.       for(i=1;i<=m;i++)
  1146.     {
  1147.       t=typ(p3[i]);if((t>5)||(t==3)) err(matqzer3);
  1148.       p2=ggcd(p2,p3[i]);
  1149.     }
  1150.       p1[j]=ldiv(p3,p2);
  1151.     }
  1152.   x=p1;
  1153.   if(gcmp0(pp))
  1154.     {
  1155.       pt=gtrans(x);
  1156.       p1=cgetg(n+1,19);
  1157.       for(j=1;j<=n;j++) p1[j]=pt[j];
  1158.       p2=det(p1);p1[n]=pt[n+1];p3=det(p1);
  1159.       p4=ggcd(p2,p3);if(!signe(p4)) err(impl,"matrixqz when the first 2 dets are zero");
  1160.       if(gcmp1(p4)) {tetpil=avma;return gerepile(av,tetpil,gcopy(x));}
  1161.       p3=factor(p4);p1=(GEN)p3[1];p2=(GEN)p3[2];nfact=lg(p1)-1;
  1162.       av1=avma;p3=cgetg(n+1,17);
  1163.       for(i=1;i<=nfact;i++)
  1164.     {
  1165.       p=(GEN)p1[i];unmodp=gmodulcp(gun,p);fl=1;
  1166.       while(fl)
  1167.         {
  1168.           pk=ker(gmul(unmodp,x));if(lg(pk)==1) fl=0;
  1169.           else
  1170.         {
  1171.           p5=centerlift(pk);p4=gdiv(gmul(x,p5),p);
  1172.           for(i=1;i<=n;i++) p3[i]=0;
  1173.           for(j=1;j<lg(p5);j++)
  1174.             {
  1175.               j1=n;while(gcmp0(coeff(p5,j1,j))) j1--;
  1176.               p3[j1]=j;
  1177.             }
  1178.           for(j=1;j<=n;j++) if(p3[j]) x[j]=p4[p3[j]];
  1179.         }
  1180.         }
  1181.       if(avma<lim) {tetpil=avma;x=gerepile(av1,tetpil,gcopy(x));}
  1182.     }
  1183.     }
  1184.   else
  1185.     {
  1186.       unmodp=gmodulcp(gun,pp);fl=1;
  1187.       while(fl)
  1188.     {
  1189.       pk=ker(gmul(unmodp,x));if(lg(pk)==1) fl=0;
  1190.       else
  1191.         {
  1192.           p5=centerlift(pk);p4=gdiv(gmul(x,p5),pp);
  1193.           for(i=1;i<=n;i++) p3[i]=0;
  1194.           for(j=1;j<lg(p5);j++)
  1195.         {
  1196.           j1=n;while(gcmp0(coeff(p5,j1,j))) j1--;
  1197.           p3[j1]=j;
  1198.         }
  1199.           for(j=1;j<=n;j++)    if(p3[j]) x[j]=p4[p3[j]];
  1200.           if(avma<lim) {tetpil=avma;x=gerepile(av1,tetpil,gcopy(x));}
  1201.         }
  1202.     }
  1203.     }
  1204.   tetpil=avma;return gerepile(av,tetpil,gcopy(x));
  1205. }
  1206.  
  1207. GEN matrixqz2(x)
  1208.      GEN x;
  1209. {
  1210.   long av=avma,tetpil,i,j,k,m,n,fl,lim,in[2];
  1211.   GEN p1;
  1212.  
  1213.   if(typ(x)!=19) err(matqzer1);
  1214.   lim=(avma+bot)>>1;
  1215.   n=lg(x)-1;if(!n) return gcopy(x);
  1216.   x=gcopy(x);
  1217.   m=lg(x[1])-1;
  1218.   for(i=1;i<=m;i++)
  1219.     {
  1220.       do
  1221.     {
  1222.       for(fl=0,j=1;(j<=n)&&(fl<2);j++)
  1223.         if(!gcmp0(coeff(x,i,j))) in[fl++]=j;
  1224.       if(fl==2)
  1225.         {
  1226.           j=(gcmp(gabs(coeff(x,i,in[0])),gabs(coeff(x,i,in[1])))>0)?in[1]:in[0];
  1227.           p1=(GEN)coeff(x,i,j);
  1228.           for(k=1;k<=n;k++)
  1229.         if(k!=j) x[k]=lsub(x[k],gmul(ground(gdiv(coeff(x,i,k),p1)),x[j]));
  1230.         }
  1231.     }
  1232.       while(fl==2);
  1233.       j=1;while((j<=n)&&gcmp0(coeff(x,i,j))) j++;
  1234.       if(j<=n) x[j]=lmul(denom(coeff(x,i,j)),x[j]);
  1235.       if(avma<lim) {tetpil=avma;x=gerepile(av,tetpil,gcopy(x));}
  1236.     }
  1237.   tetpil=avma;return gerepile(av,tetpil,hnf(x));
  1238. }
  1239.  
  1240. GEN matrixqz3(x)
  1241.      GEN x;
  1242. {
  1243.   long av=avma,av1,tetpil,i,j,j1,k,m,n,fl,lim,in[2];
  1244.   GEN p1,c,d;
  1245.  
  1246.   if(typ(x)!=19) err(matqzer1);
  1247.   lim=(avma+bot)>>1;
  1248.   n=lg(x)-1;if(!n) return gcopy(x);
  1249.   x=gcopy(x);m=lg(x[1])-1;
  1250.   c=cgeti(n+1);for(i=1;i<=n;i++) c[i]=0;
  1251.   d=cgeti(m+1);for(i=1;i<=m;i++) d[i]=0;
  1252.   av1=avma;
  1253.   for(k=1;k<=m;k++)
  1254.     {
  1255.       j=1;while((j<=n)&&(c[j]||gcmp0(coeff(x,k,j)))) j++;
  1256.       if(j<=n)
  1257.     {
  1258.       x[j]=ldiv(x[j],coeff(x,k,j));
  1259.       for(j1=1;j1<=n;j1++) if(j1!=j) x[j1]=lsub(x[j1],gmul(coeff(x,k,j1),x[j]));
  1260.       c[j]=k;d[k]=j;
  1261.     }
  1262.       if(avma<lim) {tetpil=avma;x=gerepile(av1,tetpil,gcopy(x));}
  1263.     }
  1264.   for(i=1;i<=m;i++)
  1265.     {
  1266.       do
  1267.     {
  1268.       for(fl=0,j=1;(j<=n)&&(fl<2);j++)
  1269.         if(!gcmp0(coeff(x,i,j))) in[fl++]=j;
  1270.       if(fl==2)
  1271.         {
  1272.           j=(gcmp(gabs(coeff(x,i,in[0])),gabs(coeff(x,i,in[1])))>0)?in[1]:in[0];
  1273.           p1=(GEN)coeff(x,i,j);
  1274.           for(k=1;k<=n;k++)
  1275.         if(k!=j) x[k]=lsub(x[k],gmul(ground(gdiv(coeff(x,i,k),p1)),x[j]));
  1276.         }
  1277.     }
  1278.       while(fl==2);
  1279.       j=1;while((j<=n)&&gcmp0(coeff(x,i,j))) j++;
  1280.       if(j<=n) {p1=denom(coeff(x,i,j));if(!gcmp1(p1)) x[j]=lmul(p1,x[j]);}
  1281.       if(avma<lim) {tetpil=avma;x=gerepile(av1,tetpil,gcopy(x));}
  1282.     }
  1283.   tetpil=avma;return gerepile(av,tetpil,hnf(x));
  1284. }
  1285.  
  1286. GEN kerint1(x)
  1287.      GEN x;
  1288. {
  1289.   long av=avma,tetpil;
  1290.   GEN p1,p2;
  1291.  
  1292.   p1=matrixqz3(ker(x));
  1293.   p2=lllint(p1);tetpil=avma;return gerepile(av,tetpil,gmul(p1,p2));
  1294. }
  1295.  
  1296. GEN kerint2(x)
  1297.      GEN x;
  1298. {
  1299.   long lx=lg(x), tx=typ(x),i,j,av,av1;
  1300.   GEN g,p1;
  1301.  
  1302.   if(tx!=19) err(lller1);
  1303.   av=avma;
  1304.   g=cgetg(lx,19);
  1305.   for(j=1;j<lx;j++) g[j]=lgetg(lx,18);
  1306.   for(i=1;i<lx;i++)
  1307.     for(j=1;j<=i;j++) coeff(g,i,j)=coeff(g,j,i)=(long)gscal(x[i],x[j]);
  1308.   g=lllgramall(g,1);p1=lllint(g);
  1309.   av1=avma;return gerepile(av,av1,gmul(g,p1));
  1310. }
  1311. /*
  1312. GEN oldkerint(x)
  1313.      GEN x;
  1314. {
  1315.   long lx=lg(x), tx=typ(x),i,j,av,av1;
  1316.   GEN g,p1;
  1317.  
  1318.   if(tx!=19) err(lller1);
  1319.   av=avma;
  1320.   g=cgetg(lx,19);
  1321.   for(j=1;j<lx;j++) g[j]=lgetg(lx,18);
  1322.   for(i=1;i<lx;i++)
  1323.     for(j=1;j<=i;j++) coeff(g,i,j)=coeff(g,j,i)=(long)gscal(x[i],x[j]);
  1324.   g=lllgramall0(g,1);p1=lllint(g);
  1325.   av1=avma;return gerepile(av,av1,gmul(g,p1));
  1326. }
  1327. */
  1328.  
  1329. GEN kerint(x)
  1330.      GEN x;
  1331. {
  1332.   long av=avma,av1;
  1333.   GEN g,p1;
  1334.  
  1335.   g=lllall0(x,1);p1=lllint(g);
  1336.   av1=avma;return gerepile(av,av1,gmul(g,p1));
  1337. }
  1338.  
  1339. GEN intersect(x,y)
  1340.      GEN x,y;
  1341. {
  1342.   long av=avma,tetpil,i,j,k,p;
  1343.   GEN z,p1,r;
  1344.  
  1345.   if((typ(x)!=19)||(typ(y)!=19)) err(interer1);
  1346.   if((lg(x)==1)||(lg(y)==1)) return cgetg(1,19);
  1347.   z=ker(concat(x,y));k=lg(x);p=lg(z);
  1348.   r=cgetg(p,19);
  1349.   for(j=1;j<p;j++)
  1350.     {
  1351.       p1=cgetg(k,18);r[j]=(long)p1;
  1352.       for(i=1;i<k;i++) p1[i]=coeff(z,i,j);
  1353.     }
  1354.   tetpil=avma;return gerepile(av,tetpil,gmul(x,r));
  1355. }
  1356.   
  1357.           
  1358.         
  1359.         
  1360.       
  1361.       
  1362.