home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Professional / OS2PRO194.ISO / os2 / progs / pari / pari_137 / src / alglin1.c next >
Encoding:
C/C++ Source or Header  |  1992-05-20  |  28.3 KB  |  1,132 lines
  1. /********************************************************************/
  2. /********************************************************************/
  3. /**                                                                **/
  4. /**                  ++++++++++++++++++++++++++++++                **/
  5. /**                  +                            +                **/
  6. /**                  +     ALGEBRE LINEAIRE       +                **/
  7. /**                  +                            +                **/
  8. /**                  ++++++++++++++++++++++++++++++                **/
  9. /**                                                                **/
  10. /**                        (premiere partie)                       **/
  11. /**                                                                **/
  12. /**                       copyright Babe Cool                      **/
  13. /**                                                                **/
  14. /**                                                                **/
  15. /********************************************************************/
  16. /********************************************************************/
  17.  
  18. # include "genpari.h"
  19.  
  20. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  21. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  22. /*                                                                 */
  23. /*                      TRANSPOSITION                              */
  24. /*                                                                 */
  25. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  26. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  27.  
  28. GEN     gtrans(x)
  29.      
  30.      GEN     x;
  31.      
  32. {
  33.   long    i,j,lx,tx,dx;
  34.   GEN     y,p1;
  35.   
  36.   tx=typ(x);if(tx<17) err(gtraner);
  37.   else
  38.     switch(tx)
  39.       {
  40.       case 17: y=gcopy(x);settyp(y,18);break;
  41.         
  42.       case 18: y=gcopy(x);settyp(y,17);break;
  43.         
  44.       case 19: if((lx=lg(x))==1) return cgetg(1,19);
  45.     dx=lg(x[1]);y=cgetg(dx,tx);
  46.         for(i=1;i<dx;i++)
  47.           {
  48.             p1=cgetg(lx,18);y[i]=(long)p1;
  49.             for(j=1;j<lx;j++)
  50.               p1[j]=lcopy(coeff(x,i,j));
  51.           }
  52.         break;
  53.         
  54.       default: y=gcopy(x);break;
  55.       }
  56.   return y;
  57. }
  58.  
  59. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  60. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  61. /*~                                                               ~*/
  62. /*~                  CONCATENATION ET EXTRACTION                  ~*/
  63. /*~                                                               ~*/
  64. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  65. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  66.  
  67. GEN     concat(x,y)
  68.      
  69.      GEN  x,y;
  70.      
  71. {
  72.   GEN  z,p1;
  73.   long tx=typ(x),ty=typ(y),lx=lg(x),ly=lg(y),i,dx;
  74.   
  75.   if((tx==19)&&(lx==1)) 
  76.     {
  77.       if((ty!=17)||(ly==1)) return gtomat(y);
  78.       else err(concater);
  79.     }
  80.   if((ty==19)&&(ly==1)) 
  81.     {
  82.       if((tx!=17)||(lx==1)) return gtomat(x);
  83.       else err(concater);
  84.     }
  85.   if(tx<17)
  86.     {
  87.       if(ty<17)
  88.         {
  89.           z=cgetg(3,17);z[1]=lcopy(x);
  90.           z[2]=lcopy(y);
  91.         }
  92.       else
  93.         {
  94.           if(ty!=19)
  95.             {
  96.               z=cgetg(ly+1,ty);z[1]=lcopy(x);
  97.               for(i=2;i<=ly;i++)
  98.                 z[i]=lcopy(y[i-1]);
  99.             }
  100.           else
  101.             {
  102.               if(lg(y[1])!=2) err(concater);
  103.               z=cgetg(ly+1,ty);p1=cgetg(2,18);
  104.               z[1]=(long)p1;p1[1]=lcopy(x);
  105.               for(i=2;i<=ly;i++)
  106.                 z[i]=lcopy(y[i-1]);
  107.             }
  108.         }
  109.     }
  110.   else
  111.     {
  112.       switch(tx)
  113.         {
  114.         case 17:
  115.           if(ty<17)
  116.             {
  117.               z=cgetg(lx+1,tx);z[lx]=lcopy(y);
  118.               for(i=1;i<lx;i++)
  119.                 z[i]=lcopy(x[i]);
  120.             }
  121.           else
  122.             {
  123.               switch(ty)
  124.                 {
  125.                 case 17: z=cgetg(lx+ly-1,tx);
  126.                   for(i=1;i<lx;i++)
  127.                     z[i]=lcopy(x[i]);
  128.                   for(i=1;i<ly;i++)
  129.                     z[lx+i-1]=lcopy(y[i]);
  130.                   break;
  131.                 case 18:
  132.                   if(lx==1) z=concat(x[1],y);
  133.                   else
  134.                     {
  135.                       if(ly!=1) err(concater);
  136.                       z=concat(x,y[1]);
  137.                     }
  138.                   break;
  139.                 case 19: if(lx!=ly) err(concater);
  140.                   z=cgetg(ly,ty);
  141.                   for(i=1;i<ly;i++)
  142.                     z[i]=lconcat(x[i],y[i]);
  143.                   break;
  144.                 default:;
  145.                 }
  146.             }
  147.           break;
  148.         case 18:
  149.           if(ty<17)
  150.             {
  151.               z=cgetg(lx+1,tx);z[lx]=lcopy(y);
  152.               for(i=1;i<lx;i++)
  153.                 z[i]=lcopy(x[i]);
  154.             }
  155.           else
  156.             {
  157.               switch(ty)
  158.                 {
  159.                 case 17:
  160.                   if(lx==1) z=concat(x[1],y);
  161.                   else
  162.                     {
  163.                       if(ly!=1) err(concater);
  164.                       z=concat(x,y[1]);
  165.                     }
  166.                   break;
  167.                 case 18: z=cgetg(lx+ly-1,tx);
  168.                   for(i=1;i<lx;i++)
  169.                     z[i]=lcopy(x[i]);
  170.                   for(i=1;i<ly;i++)
  171.                     z[lx+i-1]=lcopy(y[i]);
  172.                   break;
  173.                 case 19: if(lx!=lg(y[1])) err(concater);
  174.                   z=cgetg(ly+1,ty);
  175.                   z[1]=lcopy (x);
  176.                   for(i=2;i<=ly;i++)
  177.                     z[i]=lcopy(y[i-1]);
  178.                   break;
  179.                 default:;
  180.                 }
  181.             }
  182.           break;
  183.         case 19: dx=lg(x[1]);
  184.           if(ty<17)
  185.             {
  186.               if(dx!=1) err(concater);
  187.               z=cgetg(lx+1,tx);
  188.               for(i=1;i<lx;i++)
  189.                 z[i]=lcopy(x[i]);
  190.               p1=cgetg(2,18);z[lx]=(long)p1;
  191.               p1[1]=lcopy(y);
  192.             }
  193.           else
  194.             {
  195.               switch(ty)
  196.                 {
  197.                 case 17: if(lx!=ly) err(concater);
  198.                   z=cgetg(lx,tx);
  199.                   for(i=1;i<lx;i++)
  200.                     z[i]=lconcat(x[i],y[i]);
  201.                   break;
  202.                 case 18: if(dx!=ly) err(concater);
  203.                   z=cgetg(lx+1,tx);
  204.                   for(i=1;i<lx;i++)
  205.                     z[i]=lcopy(x[i]);
  206.                   z[lx]=lcopy(y);
  207.                   break;
  208.                 case 19: if(dx!=lg(y[1])) err(concater);
  209.                   z=cgetg(lx+ly-1,tx);
  210.                   for(i=1;i<lx;i++)
  211.                     z[i]=lcopy(x[i]);
  212.                   for(i=1;i<ly;i++)
  213.                     z[lx+i-1]=lcopy(y[i]);
  214.                   break;
  215.                 default:;
  216.                 }
  217.             }
  218.           break;
  219.         default:;
  220.         }
  221.     }
  222.   return z;
  223. }
  224.  
  225. GEN    vectextract(x,l)
  226.      
  227.      GEN  x,l;
  228.      
  229.      /* extraction des composantes de x suivants les bits du masque l      */
  230.      /* a usage interne donc aucune verification n'est faite. voir extract */
  231.      
  232. {
  233.   GEN  y;
  234.   long i,tx=typ(x),lx=lg(x),av,tetpil,f;
  235.   
  236.   if(!signe(l)) return cgetg(1,tx);
  237.   av=avma;i=1;
  238.   while(!mpodd(l))
  239.     {
  240.       l=shifti(l,-1);i++;
  241.     }
  242.   if(i>=lx) err(extracter3);
  243.   l=shifti(l,-1);tetpil=avma;
  244.   y=cgetg(2,tx);y[1]=lcopy(x[i]);
  245.   i++;
  246.   while(!gcmp0(l)&&(i<lx))
  247.     {
  248.       f=mpodd(l);l=shifti(l,-1);tetpil=avma;
  249.       if(f) y=concat(y,x[i]);
  250.       i++;
  251.     }
  252.   if(!gcmp0(l)) err(extracter3);
  253.   y=gerepile(av,tetpil,y);
  254.   return y;
  255. }
  256.  
  257. GEN    extract(x,l)
  258.      
  259.      GEN   x,l;
  260.      
  261. {
  262.   long tl=typ(l),ll,lx,i,tx=typ(x),in;
  263.   GEN y;
  264.   
  265.   if(tx<17) err(extracter1);
  266.   if(tl==1) return vectextract(x,l);
  267.   if((tl==17)||(tl==18))
  268.     {
  269.       ll=lg(l);y=cgetg(ll,tx);lx=lg(x);
  270.       for(i=1;i<ll;i++)
  271.         {
  272.           in=itos(l[i]);if((in>=lx)||(in<=0)) err(extracter3);
  273.           y[i]=lcopy(x[in]);
  274.         }
  275.       return y;
  276.     }
  277.   err(extracter2);
  278. }
  279.  
  280. GEN    matextract(x,l1,l2)
  281.      
  282.      GEN   x,l1,l2;
  283.      
  284. {
  285.   GEN  y;
  286.   long av,tetpil;
  287.   
  288.   if(typ(x)!=19) err(matextrer);
  289.   av=avma;y=extract(gtrans(extract(x,l2)),l1);tetpil=avma;
  290.   return gerepile(av,tetpil,gtrans(y));
  291. }
  292.  
  293.  
  294. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  295. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  296. /*                                                                 */
  297. /*              OPERATIONS SCALAIRES-MATRICES                      */
  298. /*                                                                 */
  299. /*                        ET DIVERS                                */
  300. /*                                                                 */
  301. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  302. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  303.  
  304.  
  305. GEN     gscalmat(x,n)       /* cree la matrice carree n X n */
  306.      /* contenant x*I                */
  307.      GEN     x;
  308.      long    n;
  309.      
  310. {
  311.   long    i,j,z;
  312.   GEN     y;
  313.   
  314.   z = lcopy(x);
  315.   y=cgetg(n+1,19);
  316.   for(i=1;i<=n;i++)
  317.     {
  318.       y[i]=lgetg(n+1,18);
  319.       for(j=1;j<=n;j++)
  320.         coeff(y,j,i)=(i==j ? z : zero);
  321.     }
  322.   return y;
  323. }
  324.  
  325. GEN     gscalsmat(x,n)      /* idem au precedent avec x long du C   */
  326.      
  327.      long    x,n;
  328.      
  329. {
  330.   long    i,j,z;
  331.   GEN     y;
  332.   
  333.   z=lstoi(x);
  334.   y=cgetg(n+1,19);
  335.   for(i=1;i<=n;i++)
  336.     {
  337.       y[i]=lgetg(n+1,18);
  338.       for(j=1;j<=n;j++)
  339.         coeff(y,j,i)=(i==j ? z : zero);
  340.     }
  341.   return y;
  342. }
  343.  
  344. GEN   idmat(n)
  345.      
  346.      long  n;
  347. {
  348.   return gscalmat(gun,n);
  349. }
  350.  
  351. GEN   gtomat(x)
  352.      
  353.      GEN  x;
  354.      
  355. {
  356.   GEN  y,p1;
  357.   long tx=typ(x),lx,i;
  358.   
  359.   if(tx<17)
  360.     {
  361.       y=cgetg(2,19);p1=cgetg(2,18);y[1]=(long)p1;
  362.       p1[1]=lcopy(x);
  363.     }
  364.   else
  365.     switch(tx)
  366.       {
  367.       case 17: lx=lg(x);y=cgetg(lx,19);
  368.         for(i=1;i<lx;i++)
  369.           {
  370.             p1=cgetg(2,18);p1[1]=lcopy(x[i]);
  371.             y[i]=(long)p1;
  372.           }
  373.         break;
  374.       case 18: y=cgetg(2,19);y[1]=lcopy(x);break;
  375.       case 19: y=gcopy(x);break;
  376.       }
  377.   return y;
  378. }
  379.  
  380. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  381. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  382. /*                                                                 */
  383. /*                   ADDITION SCALAIRE +  MATRICE                  */
  384. /*                                                                 */
  385. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  386. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  387.  
  388. GEN     gaddmat(x,y)        /* cree la matrice carree contenant x*I+y     */
  389.      
  390.      GEN     x,y;
  391.      
  392. {
  393.   long    ly,dy,i,j;
  394.   GEN     z;
  395.   
  396.   ly=lg(y);dy=lg(y[1]);
  397.   if((typ(y)!=19) || (ly!=dy)) err(gadmaer);
  398.   z=cgetg(ly,19);
  399.   for(i=1;i<ly;i++)
  400.     {
  401.       z[i]=lgetg(dy,18);
  402.       for(j=1;j<dy;j++)
  403.     coeff(z,j,i)=(i==j ? ladd(x,coeff(y,j,i)) : lcopy(coeff(y,j,i)));
  404.     }
  405.   return z;
  406. }
  407.  
  408. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  409. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  410. /*                                                                 */
  411. /*                      ADDITION SHORT +  MATRICE                  */
  412. /*                                                                 */
  413. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  414. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  415.  
  416. GEN     gaddsmat(s,y)       /* idem au precedent avec x long du C   */
  417.      
  418.      long    s;
  419.      GEN     y;
  420.      
  421. {
  422.   long    ly,dy,i,j;
  423.   GEN     z;
  424.   
  425.   ly=lg(y);dy=lg(y[1]);
  426.   if((typ(y)!=19) || (ly!=dy)) err(gadsmaer);
  427.   z=cgetg(ly,19);
  428.   for(i=1;i<ly;i++)
  429.     {
  430.       z[i]=lgetg(dy,18);
  431.       for(j=1;j<dy;j++)
  432.     coeff(z,j,i)=(i==j ? laddsg(s,coeff(y,j,i)) : lcopy(coeff(y,j,i)));
  433.     }
  434.   return z;
  435. }
  436.  
  437. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  438. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  439. /*                                                                 */
  440. /*                      RESOLUTION DE A X=B                        */
  441. /*                      (METHODE DE GAUSS)                         */
  442. /*                                                                 */
  443. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  444. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  445.  
  446. GEN     gauss(a,b)
  447.      
  448.      GEN     a,b;
  449.      
  450. {
  451.   long  p,u,m,nbli,nbco,i,j,k,av1,av2,av3,av4;
  452.   GEN   aa,x;
  453.   
  454.   if(typ(b)==19) return invmulmat(a,b);
  455.   nbco=lg(a)-1;nbli=lg(a[1])-1;
  456.   if (nbco!=nbli) err(gausser1);
  457.   x=cgetg(nbli+1,18);av1=avma;
  458.   for (j=1;j<=nbco;j++) x[j]=b[j];
  459.   aa=cgetg(nbco+1,19);
  460.   for (j=1;j<=nbco;j++)
  461.     {
  462.       aa[j]=lgetg(nbli+1,18);
  463.       for (i=1;i<=nbli;i++) coeff(aa,i,j)=coeff(a,i,j);
  464.     }
  465.   for (i=1;i<nbli;i++)
  466.     {
  467.       p=coeff(aa,i,i);k=i;
  468.       if (gcmp0(p))
  469.     {
  470.       for (k=i+1;(k<=nbli)&&gcmp0(coeff(aa,k,i));k++);
  471.       if (k>nbco) err(matinv1);
  472.       else
  473.         {
  474.           for (j=i;j<=nbco;j++)
  475.         {
  476.           u=coeff(aa,i,j);coeff(aa,i,j)=coeff(aa,k,j);
  477.           coeff(aa,k,j)=u;
  478.         }
  479.           u=x[i];x[i]=x[k];x[k]=u;
  480.           p=coeff(aa,i,i);
  481.         }
  482.     }
  483.       for (k=i+1;k<=nbli;k++)
  484.     {
  485.       m=coeff(aa,k,i);
  486.       if (!gcmp0(m))
  487.         {
  488.           m=ldiv(m,p);
  489.           for (j=i+1;j<=nbco;j++)
  490.         coeff(aa,k,j)=lsub(coeff(aa,k,j),gmul(m ,coeff(aa,i,j)));
  491.           x[k]=lsub(x[k],gmul(m,x[i]));
  492.         }
  493.     }
  494.     }
  495.   
  496.   /* Resolution systeme triangularise */
  497.   av2=avma;
  498.   p=coeff(aa,nbli,nbco);
  499.   if (gcmp0(p)) err(matinv1);
  500.   else
  501.     {
  502.       x[nbli]=ldiv(x[nbli],p);
  503.       for (i=nbli-1;i>0;i--)
  504.     {
  505.       av3=avma;m=x[i];
  506.       for (j=i+1;j<=nbco;j++)
  507.         m= lsub(m,gmul(coeff(aa,i,j),x[j]));
  508.       av4=avma;
  509.       x[i]=lpile(av3,av4,gdiv(m,coeff(aa,i,i)));
  510.     }
  511.     }
  512.   gerepile(av1,av2,1);
  513.   return x;
  514. }
  515.  
  516. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  517. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  518. /*                                                                 */
  519. /*            RANG D'UNE MATRICE m lignes x n colonnes             */
  520. /*                                                                 */
  521. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  522. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  523.  
  524. long rank(x)
  525.      GEN  x;
  526.      
  527. {
  528.   GEN c,mun,p;
  529.   long i,j,k,r,t,n,m,av;
  530.   
  531.   if (typ(x)!=19) err(kerer1);
  532.   r=n=lg(x)-1;if(!r) return r;
  533.   m=lg(x[1])-1;av=avma;
  534.   x=gcopy(x);c=cgeti(m+1);
  535.   mun=gneg(un);
  536.   for(k=1;k<=m;k++) c[k]=0;
  537.   for(k=1;k<=n;k++)
  538.     {
  539.       j=1;
  540.       while((j<=m)&&(c[j]||gcmp0(coeff(x,j,k)))) j++;
  541.       if (j<=m)
  542.     {
  543.       p=gdivsg(-1,coeff(x,j,k));
  544.       coeff(x,j,k)=(long)mun;
  545.       for(i=k+1;i<=n;i++) coeff(x,j,i)=lmul(p,coeff(x,j,i));
  546.       for(t=1;t<=m;t++)
  547.         if(t!=j)
  548.           {
  549.         p=(GEN)coeff(x,t,k);
  550.         for(i=k+1;i<=n;i++) coeff(x,t,i)=ladd(coeff(x,t,i),gmul(p,coeff(x,j,i)));
  551.         coeff(x,t,k)=zero;
  552.           }
  553.       c[j]=k;
  554.     }                  
  555.       else r--;
  556.     }
  557.   avma=av;
  558.   return r;
  559. }
  560.  
  561. GEN indexrank(x)
  562.      GEN  x;
  563.      
  564. {
  565.   GEN c,d,mun,p,y,p1,p2;
  566.   long i,j,k,r,t,n,m,av,tetpil;
  567.   
  568.   if (typ(x)!=19) err(kerer1);
  569.   r=n=lg(x)-1;if(!r) {y=cgetg(3,17);y[1]=lgetg(1,17);y[2]=lgetg(1,17);return y;}
  570.   m=lg(x[1])-1;av=avma;
  571.   x=gcopy(x);c=cgeti(m+1);d=cgeti(n+1);
  572.   mun=gneg(un);
  573.   for(k=1;k<=m;k++) c[k]=0;
  574.   for(k=1;k<=n;k++)
  575.     {
  576.       j=1;
  577.       while((j<=m)&&(c[j]||gcmp0(coeff(x,j,k)))) j++;
  578.       if (j<=m)
  579.     {
  580.       p=gdivsg(-1,coeff(x,j,k));
  581.       coeff(x,j,k)=(long)mun;
  582.       for(i=k+1;i<=n;i++) coeff(x,j,i)=lmul(p,coeff(x,j,i));
  583.       for(t=1;t<=m;t++)
  584.         if(t!=j)
  585.           {
  586.         p=(GEN)coeff(x,t,k);
  587.         for(i=k+1;i<=n;i++) coeff(x,t,i)=ladd(coeff(x,t,i),gmul(p,coeff(x,j,i)));
  588.         coeff(x,t,k)=zero;
  589.           }
  590.       c[j]=k;d[k]=j;
  591.     }                  
  592.       else {r--;d[k]=0;}
  593.     }
  594.   p1=cgetg(r+1,17);p2=cgetg(r+1,17);
  595.   for(i=0,k=1;k<=n;k++) if(d[k]) {p1[++i]=lstoi(d[k]);p2[i]=lstoi(k);}
  596.   tetpil=avma;y=cgetg(3,17);y[1]=(long)sort(p1);y[2]=lcopy(p2);
  597.   return gerepile(av,tetpil,y);
  598. }
  599.  
  600. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  601. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  602. /*                                                                 */
  603. /*          NOYAU D'UNE MATRICE m lignes x n colonnes              */
  604. /*   ( Retourne une matrice de n-rang vecteurs independants )      */
  605. /*                                                                 */
  606. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  607. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  608.  
  609. GEN     keri(x)     /* Programme pour types ENTIERS */
  610.      GEN  x;
  611.      
  612. {
  613.   GEN c,d,y,v,pp;
  614.   long i,j,k,r,t,n,n1,m,nbmot,av,av1;
  615.   long p,p0,p1,q;
  616.   
  617.   if (typ(x)!=19) err(kerer1);
  618.   n1=lg(x);n=n1-1;if(!n) return cgetg(1,19);
  619.   m=lg(x[1])-1;av=avma;
  620.   nbmot=200;
  621.   c=cgetg(n1,19);
  622.   for(j=1;j<=n;j++)
  623.     {
  624.       p=c[j]=lgetg(m+1,18);
  625.       for(i=1;i<=m;i++)
  626.     affii(coeff(x,i,j),((GEN)p)[i]=lgeti(nbmot));
  627.     }
  628.   x=c;
  629.   p=un;
  630.   pp=cgetg(n+1,18);
  631.   for(j=1;j<=n;j++) pp[j]=lgeti(nbmot);
  632.   c=cgeti(m+1);for(k=1;k<=m;k++) c[k]=0;
  633.   d=cgeti(n1);
  634.   av1=avma;
  635.   for(r=0,k=1;k<=n;k++)
  636.     {
  637.       j=1;
  638.       while((j<=m)&&(c[j]||!signe(coeff(x,j,k)))) j++;
  639.       if (j<=m)
  640.     {
  641.       p0=p;p=coeff(x,j,k);
  642.       for(t=1;t<=m;t++)
  643.         if(t!=j)
  644.           {
  645.         q=coeff(x,t,k);
  646.         for(i=k+1;i<=n;i++)
  647.           {
  648.             p1=lsubii(mulii(p,coeff(x,t,i)),mulii(q,coeff(x,j,i)));
  649.             if(k>1) diviiz(p1,p0,coeff(x,t,i));
  650.             else affii(p1,coeff(x,t,i));
  651.           }
  652.           }
  653.       c[j]=k;d[k]=j;
  654.       avma=av1;
  655.     }          
  656.        else {r++;d[k]=0;affii(p,pp[k]);}
  657.  
  658.     }
  659.   if(r)    /* Il y a un noyau non nul */
  660.     {
  661.       av1=avma;
  662.       y=cgetg(r+1,19);
  663.       for(j=k=1;j<=r;j++,k++)
  664.     {
  665.       while(d[k]) k++;
  666.       y[j]=(long)(v=cgetg(n1,18));
  667.       for(i=1;i<k;i++) v[i]= d[i]? lcopy(coeff(x,d[i],k)) : zero;
  668.       v[k]=lnegi(pp[k]);
  669.       for(i=k+1;i<=n;i++) v[i]=zero;
  670.     }
  671.       return gerepile(av,av1,y);
  672.     }
  673.   else {avma=av;y=cgetg(1,19);return y;}
  674. }
  675.  
  676. GEN     ker(x)     /* Programme pour types exacts */
  677.      GEN  x;
  678.      
  679. {
  680.   GEN c,d,y,mun,p;
  681.   long i,j,k,r,t,n,n1,m,av,av1,av2;
  682.   
  683.   if (typ(x)!=19) err(kerer1);
  684.   n1=lg(x);n=n1-1;if(!n) return cgetg(1,19);
  685.   m=lg(x[1])-1;av=avma;x=gcopy(x);mun=gneg(un);r=0;
  686.   c=cgeti(m+1);for(k=1;k<=m;k++) c[k]=0;
  687.   d=cgeti(n1);
  688.   av1=avma;
  689.   for(k=1;k<=n;k++)
  690.     {
  691.       j=1;
  692.       while((j<=m)&&(c[j]||gcmp0(coeff(x,j,k)))) j++;
  693.       if (j<=m)
  694.     {
  695.       
  696.       p=gdivsg(-1,coeff(x,j,k));
  697.       coeff(x,j,k)=(long)mun;
  698.       for(i=k+1;i<=n;i++) coeff(x,j,i)=lmul(p,coeff(x,j,i));
  699.       for(t=1;t<=m;t++)
  700.         if(t!=j)
  701.           {
  702.         p=(GEN)coeff(x,t,k);
  703.         for(i=k+1;i<=n;i++) coeff(x,t,i)=ladd(coeff(x,t,i),gmul(p,coeff(x,j,i)));
  704.         coeff(x,t,k)=zero;
  705.           }
  706.       c[j]=k;d[k]=j;
  707.       av2=avma;
  708.       x=gerepile(av1,av2,gcopy(x));
  709.     }          
  710.       else {r++;d[k]=0;}
  711.     }
  712.   if(r)
  713.     {
  714.       av1=avma;
  715.       y=cgetg(r+1,19);
  716.       for(j=k=1;j<=r;j++,k++)
  717.     {
  718.       while(d[k]) k++;
  719.       y[j]=(long)(p=cgetg(n1,18));
  720.       for(i=1;i<k;i++) p[i]=d[i]? lcopy(coeff(x,d[i],k)):zero;
  721.       p[k]=un;
  722.       for(i=k+1;i<=n;i++) p[i]=zero;
  723.     }
  724.       return gerepile(av,av1,y);
  725.     }
  726.   else {avma=av;y=cgetg(1,19);return y;}
  727. }
  728.  
  729. GEN     image(x)     /* Programme pour types exacts */
  730.      GEN  x;
  731.      
  732. {
  733.   GEN c,d,y,mun,p,x1;
  734.   long i,j,k,r,t,n,n1,m,av,av1,av2;
  735.   
  736.   if (typ(x)!=19) err(kerer1);
  737.   n1=lg(x);n=n1-1;if(!n) return cgetg(1,19);
  738.   m=lg(x[1])-1;av=avma;x1=gcopy(x);mun=gneg(un);r=0;
  739.   c=cgeti(m+1);for(k=1;k<=m;k++) c[k]=0;
  740.   d=cgeti(n1);
  741.   av1=avma;
  742.   for(k=1;k<=n;k++)
  743.     {
  744.       j=1;
  745.       while((j<=m)&&(c[j]||gcmp0(coeff(x1,j,k)))) j++;
  746.       if (j<=m)
  747.     {
  748.       
  749.       p=gdivsg(-1,coeff(x1,j,k));
  750.       coeff(x1,j,k)=(long)mun;
  751.       for(i=k+1;i<=n;i++) coeff(x1,j,i)=lmul(p,coeff(x1,j,i));
  752.       for(t=1;t<=m;t++)
  753.         if(t!=j)
  754.           {
  755.         p=(GEN)coeff(x1,t,k);
  756.         for(i=k+1;i<=n;i++) coeff(x1,t,i)=ladd(coeff(x1,t,i),gmul(p,coeff(x1,j,i)));
  757.         coeff(x1,t,k)=zero;
  758.           }
  759.       c[j]=k;d[k]=j;
  760.       av2=avma;
  761.       x1=gerepile(av1,av2,gcopy(x1));
  762.     }          
  763.       else {r++;d[k]=0;}
  764.     }
  765.   if(r)
  766.     {
  767.       av1=avma;
  768.       y=cgetg(n-r+1,19);
  769.       for(j=k=1;j<=n-r;j++,k++)
  770.     {
  771.       while(!d[k]) k++;
  772.       y[j]=lcopy(x[k]);
  773.     }
  774.       return gerepile(av,av1,y);
  775.     }
  776.   else {avma=av;return gcopy(x);}
  777. }
  778.  
  779. GEN inverseimage(mat,y)
  780.      GEN mat,y;
  781. {
  782.   long av=avma,nbcol,i,j,l,tetpil;
  783.   GEN met,noyau,lastcoeff,invimag;
  784.  
  785.   if ((typ(mat)!=19)||(typ(y)!=18)) err(kerer1);
  786.   nbcol=lg(mat);
  787.   met=cgetg(nbcol+1,19);
  788.   for(j=1;j<=nbcol-1;j++) met[j]=mat[j];
  789.   met[nbcol]=(long)y;
  790.   noyau=ker(met);l=lg(noyau)-1;
  791.   if(!l) {avma=av;return cgetg(1,18);}
  792.   lastcoeff=gneg(coeff(noyau,nbcol,l));
  793.   if(gcmp0(lastcoeff)) {avma=av;return cgetg(1,18);}
  794.   tetpil=avma;
  795.   invimag=cgetg(nbcol,18);
  796.   for(i=1;i<=nbcol-1;i++) invimag[i]=ldiv(coeff(noyau,i,l),lastcoeff);
  797.   return gerepile(av,tetpil,invimag);
  798. }
  799.  
  800. GEN     kerreel(x,prec)    /* Programme pour types non exacts    */
  801.      /* gestion de pile a la fin seulement */
  802.      GEN  x;
  803.      long prec;
  804. {
  805.   GEN c,d,y,mun,p,eps;
  806.   long i,j,k,r,t,n,n1,m,av,av1;
  807.   
  808.   if (typ(x)!=19) err(kerer1);
  809.   n1=lg(x);n=n1-1;if(!n) return cgetg(1,19);
  810.   m=lg(x[1])-1;av=avma;
  811.   eps=cgetr(3);eps[2]=0x80000000;eps[1]=0x01800010-((prec-2)<<5);
  812.   x=gcopy(x);
  813.   mun=gneg(un);r=0;
  814.   c=cgeti(m+1);for(k=1;k<=m;k++) c[k]=0;
  815.   d=cgeti(n1);
  816.   for(k=1;k<=n;k++)
  817.     {
  818.       j=1;
  819.       while((j<=m)&&(c[j]||(gcmp(gabs(coeff(x,j,k),5),eps)<0))) j++;
  820.       if (j<=m)
  821.     {
  822.       
  823.       p=gdivsg(-1,coeff(x,j,k));
  824.       coeff(x,j,k)=(long)mun;
  825.       for(i=k+1;i<=n;i++) coeff(x,j,i)=lmul(p,coeff(x,j,i));
  826.       for(t=1;t<=m;t++)
  827.         if(t!=j)
  828.           {
  829.         p=(GEN)coeff(x,t,k);
  830.         for(i=k+1;i<=n;i++) coeff(x,t,i)=ladd(coeff(x,t,i),gmul(p,coeff(x,j,i)));
  831.         coeff(x,t,k)=zero;
  832.           }
  833.       c[j]=k;d[k]=j;
  834.     }                  
  835.       else{r++;d[k]=0;}
  836.     }
  837.   if(r)
  838.     {
  839.       av1=avma;
  840.       y=cgetg(r+1,19);
  841.       for(j=k=1;j<=r;j++,k++)
  842.     {
  843.       while(d[k]) k++;
  844.       y[j]=(long)(p=cgetg(n1,18));
  845.       for(i=1;i<k;i++) p[i]=d[i]? lcopy(coeff(x,d[i],k)):zero;
  846.       p[k]=un;
  847.       for(i=k+1;i<=n;i++) p[i]=zero;
  848.     }
  849.       return gerepile(av,av1,y);
  850.     }
  851.   else {avma=av;y=cgetg(1,19);return y;}
  852. }
  853.  
  854. /* Etant donnee une matrice nxk de rang k<=n, on trouve une matrice nxn
  855. inversible dont les k premieres colonnes forment la matrice initiale;
  856. on ne verifie pas que les k colonnes sont lineairement independantes. */
  857.  
  858. GEN suppl(x)
  859.      GEN x;
  860. {
  861.   long av=avma,tetpil,k,n,s,t;
  862.   GEN y,p1,p2;
  863.  
  864.   if(typ(x)!=19) err(kerer1);
  865.   k=lg(x)-1;if(!k) err(suppler1);
  866.   n=lg(x[1])-1;if(k>n) err(suppler2);
  867.   s=0;y=idmat(n);
  868.   while(s<k)
  869.     {
  870.       s++;p1=gauss(y,x[s]);t=s;
  871.       while((t<=n)&&gcmp0(p1[t])) t++;
  872.       if(t>n) err(suppler2);
  873.       p2=(GEN)y[s];y[s]=x[s];if(s!=t) y[t]=(long)p2;
  874.     }
  875.   tetpil=avma;return gerepile(av,tetpil,gcopy(y));
  876. }
  877.  
  878. GEN image2(x)
  879.      GEN x;
  880. {
  881.   long av=avma,tetpil,k,n,i;
  882.   GEN p1,p2;
  883.  
  884.   if(typ(x)!=19) err(kerer1);
  885.   k=lg(x)-1;if(!k) return gcopy(x);
  886.   n=lg(x[1])-1;p1=ker(x);k=lg(p1)-1;
  887.   if(k) p1=suppl(p1);else p1=idmat(n);
  888.   n=lg(p1)-1;
  889.   tetpil=avma;p2=cgetg(n-k+1,19);
  890.   for(i=k+1;i<=n;i++) p2[i-k]=lmul(x,p1[i]);
  891.   return gerepile(av,tetpil,p2);
  892. }
  893.  
  894. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  895. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  896. /*                                                                 */
  897. /*                         VECTEURS PROPRES                        */
  898. /*            (matrice de vecteurs propres independants            */
  899. /*             classes par valeurs propres croissantes )           */
  900. /*                                                                 */
  901. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  902. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  903.  
  904. GEN eigen(x,prec)
  905.      GEN x;
  906.      long prec;
  907.      
  908. {
  909.   GEN y,z,rr,p,ssesp,eps;
  910.   long j,k,n,ly,av,av1,nbrac,nk,r1,r2,r3,flag;
  911.   
  912.   
  913.   n=lg(x);
  914.   av=avma;
  915.   eps=cgetr(3);eps[2]=0x80000000;eps[1]=0x01800010-((prec-2)<<5);
  916.   y=cgetg(n,19);ly=1;
  917.   z=gcopy(x);
  918.   p=caradj(x,0,0);rr=roots(p,prec);nbrac=lg(rr)-1;
  919. /* Bien sur ce n'est pas comme cela qu'on doit calculer les valeurs propres !*/
  920.   for(k=1;k<=nbrac;k++)
  921.     {
  922.       r2=rr[k];flag=0;
  923.       if(k>1) if(gcmp(gabs(gsub(r1,r2),5),eps)>0) flag=1;
  924.       if(flag||(k==1))
  925.     {
  926.       r3=lround(r2);if(gcmp(gabs(gsub(r2,r3),5),eps)<0) r2=r3;
  927.       {
  928.         for(j=1;j<n;j++) coeff(z,j,j)=lsub(coeff(x,j,j),r2);
  929.         ssesp=kerreel(z,prec);
  930.         nk=lg(ssesp)-1;
  931.         for(j=1;j<=nk;j++,ly++) y[ly]=ssesp[j];
  932.       }
  933.       r1=r2;
  934.     }
  935.     }
  936.   z=cgetg(ly,19);
  937.   av1=avma;
  938.   for(k=1;k<ly;k++) z[k]=y[k];
  939.   return gerepile(av,av1,gcopy(z));
  940. }
  941. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  942. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  943. /*                                                                 */
  944. /*                           DETERMINANT                           */
  945. /*                                                                 */
  946. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  947. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  948.  
  949. /* ===================================================================*/
  950. /*       Determinant types exacts : 1er pivot non nul                 */
  951. /*--------------------------------------------------------------------*/
  952.  
  953. GEN     det2(a)
  954.      
  955.      GEN     a;
  956.      
  957. {
  958.   long  p,u,m,nbli,nbco,i,j,k,av,av1,s;
  959.   GEN   aa,p1,x;
  960.   
  961.   if (typ(a)!=19) err(mattype1);
  962.   nbco=lg(a)-1;nbli=lg(a[1])-1;
  963.   if (nbco!=nbli) err(mattype1);
  964.   av=avma;x=gun;s=1;
  965.   aa=cgetg(nbco+1,19);
  966.   
  967.   for (j=1;j<=nbco;j++)
  968.     {
  969.       aa[j]=lgetg(nbli+1,18);
  970.       for (i=1;i<=nbli;i++) coeff(aa,i,j)=coeff(a,i,j);
  971.     }
  972.   
  973.   for (i=1;i<nbco;i++)
  974.     {
  975.       p=coeff(aa,i,i);k=i;
  976.       if(gcmp0(p))
  977.     {
  978.       
  979.       for (k=i+1;(k<=nbco)&&gcmp0(coeff(aa,i,k));k++);
  980.       if (k>nbco)
  981.         {
  982.           avma=av;return gzero;
  983.         }
  984.       else
  985.         {
  986.           p=coeff(aa,i,k);
  987.           u=aa[k];aa[k]=aa[i];aa[i]=u;
  988.           s= -s;
  989.         }
  990.     }
  991.       x=gmul(x,p);
  992.       
  993.       for (k=i+1;k<=nbco;k++)
  994.     {
  995.       m=coeff(aa,i,k);
  996.       if (!gcmp0(m))
  997.         {
  998.           m=ldiv(m,p);
  999.           for (j=i+1;j<=nbli;j++)
  1000.         {
  1001.           p1=gmul(m,coeff(aa,j,i));
  1002.           coeff(aa,j,k)=lsub(coeff(aa,j,k),p1);
  1003.         }
  1004.         }
  1005.     }
  1006.     }
  1007.   if(s<0) x=gneg(x);
  1008.   av1=avma;
  1009.   return gerepile(av,av1,gmul(x,coeff(aa,nbli,nbco)));
  1010. }
  1011.  
  1012. /* ===================================================================*/
  1013. /*     Determinant dans un anneau A : Tous les calculs dans A         */
  1014. /*     division par le pivot precedent ( methode de Bareiss)          */
  1015. /*--------------------------------------------------------------------*/
  1016.  
  1017. GEN     det(a)
  1018.      
  1019.      GEN     a;
  1020.      
  1021. {
  1022.   long  p,pprec,u,m,nbli,nbco,i,j,k,av,av1,s;
  1023.   GEN   aa,p1;
  1024.   
  1025.   
  1026.   if (typ(a)!=19) err(mattype1);
  1027.   nbco=lg(a)-1;nbli=lg(a[1])-1;
  1028.   if (nbco!=nbli) err(mattype1);
  1029.   av=avma;
  1030.   aa=cgetg(nbco+1,19);
  1031.   
  1032.   for (j=1;j<=nbco;j++)
  1033.     {
  1034.       aa[j]=lgetg(nbli+1,18);
  1035.       for (i=1;i<=nbli;i++) coeff(aa,i,j)=coeff(a,i,j);
  1036.     }
  1037.   pprec=un;s=1;
  1038.   for (i=1;i<nbco;i++)
  1039.     {
  1040.       p=coeff(aa,i,i);k=i;
  1041.       if(gcmp0(p))
  1042.         {
  1043.           
  1044.           for (k=i+1;(k<=nbco)&&gcmp0(coeff(aa,i,k));k++);
  1045.           if (k>nbco)
  1046.             {
  1047.               avma=av;return gzero;
  1048.             }
  1049.           else
  1050.             {
  1051.               p=coeff(aa,i,k);
  1052.               u=aa[k];aa[k]=aa[i];aa[i]=u;
  1053.               s= -s;
  1054.             }
  1055.         }
  1056.       for (k=i+1;k<=nbco;k++)
  1057.         {
  1058.           m=coeff(aa,i,k);
  1059.           for (j=i+1;j<=nbli;j++)
  1060.             {
  1061.               p1=gsub(gmul(p,coeff(aa,j,k)),gmul(m,coeff(aa,j,i)));
  1062.           if((typ(p1)==10)&&(typ(pprec)==10)&&(varn(p1)==varn(pprec)))
  1063.         coeff(aa,j,k)=ldeuc(p1,pprec);
  1064.               else coeff(aa,j,k)=ldiv(p1,pprec);
  1065.             }
  1066.         }
  1067.       pprec=p;
  1068.     }
  1069.   av1=avma;
  1070.   return (s>0) ? gerepile(av,av1,gcopy(coeff(aa,nbli,nbco))) : gerepile(av,av1,gneg(coeff(aa,nbli,nbco)));
  1071. }
  1072.  
  1073. /* ===================================================================*/
  1074. /*              Determinant reel : pivot maximal                      */
  1075. /*--------------------------------------------------------------------*/
  1076.  
  1077.  
  1078. GEN     detreel(a)
  1079.      
  1080.      GEN     a;
  1081.      
  1082. {
  1083.   long  p,u,m,nbli,nbco,i,j,k,av,av1,s;
  1084.   GEN   aa,p1,x;
  1085.   
  1086.   if (typ(a)!=19) err(mattype1);
  1087.   nbco=lg(a)-1;nbli=lg(a[1])-1;
  1088.   if (nbco!=nbli) err(mattype1);
  1089.   av=avma;s=1;x=gun;
  1090.   aa=cgetg(nbco+1,19);
  1091.   
  1092.   for (j=1;j<=nbco;j++)
  1093.     {
  1094.       aa[j]=lgetg(nbli+1,18);
  1095.       for (i=1;i<=nbli;i++) coeff(aa,i,j)=coeff(a,i,j);
  1096.     }
  1097.   
  1098.   for (i=1;i<nbco;i++)
  1099.     {
  1100.       p=labs(coeff(aa,i,i));k=i;
  1101.       for(j=i+1;j<=nbco;j++)
  1102.     if(gcmp(p1=gabs(coeff(aa,i,j),3),p)>0) {p=(long)p1;k=j;}
  1103.       if(gcmp0(p))
  1104.     {
  1105.       av1=av;return gerepile(av,av1,gcopy(p));
  1106.     }
  1107.       else
  1108.     {
  1109.       p=coeff(aa,i,k);
  1110.       if(k>i)
  1111.         {
  1112.           u=aa[k];aa[k]=aa[i];aa[i]=u;
  1113.           s= -s;
  1114.         }
  1115.     }
  1116.       x=gmul(x,p);
  1117.       
  1118.       for (k=i+1;k<=nbco;k++)
  1119.     {
  1120.       m=coeff(aa,i,k);
  1121.       if (!gcmp0(m))
  1122.         {
  1123.           m=ldiv(m,p);
  1124.           for (j=i+1;j<=nbli;j++)
  1125.         coeff(aa,j,k)=lsub(coeff(aa,j,k),gmul(m,coeff(aa,j,i)));
  1126.         }
  1127.     }
  1128.     }
  1129.   if(s<0) x=gneg(x);
  1130.   av1=avma;return gerepile(av,av1,gmul(x,coeff(aa,nbli,nbco)));
  1131. }
  1132.