home *** CD-ROM | disk | FTP | other *** search
- /********************************************************************/
- /********************************************************************/
- /** **/
- /** ++++++++++++++++++++++++++++++ **/
- /** + + **/
- /** + ALGEBRE LINEAIRE + **/
- /** + + **/
- /** ++++++++++++++++++++++++++++++ **/
- /** **/
- /** (premiere partie) **/
- /** **/
- /** copyright Babe Cool **/
- /** **/
- /** **/
- /********************************************************************/
- /********************************************************************/
-
- # include "genpari.h"
-
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /* */
- /* TRANSPOSITION */
- /* */
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
-
- GEN gtrans(x)
-
- GEN x;
-
- {
- long i,j,lx,tx,dx;
- GEN y,p1;
-
- tx=typ(x);if(tx<17) err(gtraner);
- else
- switch(tx)
- {
- case 17: y=gcopy(x);settyp(y,18);break;
-
- case 18: y=gcopy(x);settyp(y,17);break;
-
- case 19: if((lx=lg(x))==1) return cgetg(1,19);
- dx=lg(x[1]);y=cgetg(dx,tx);
- for(i=1;i<dx;i++)
- {
- p1=cgetg(lx,18);y[i]=(long)p1;
- for(j=1;j<lx;j++)
- p1[j]=lcopy(coeff(x,i,j));
- }
- break;
-
- default: y=gcopy(x);break;
- }
- return y;
- }
-
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~ ~*/
- /*~ CONCATENATION ET EXTRACTION ~*/
- /*~ ~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
-
- GEN concat(x,y)
-
- GEN x,y;
-
- {
- GEN z,p1;
- long tx=typ(x),ty=typ(y),lx=lg(x),ly=lg(y),i,dx;
-
- if((tx==19)&&(lx==1))
- {
- if((ty!=17)||(ly==1)) return gtomat(y);
- else err(concater);
- }
- if((ty==19)&&(ly==1))
- {
- if((tx!=17)||(lx==1)) return gtomat(x);
- else err(concater);
- }
- if(tx<17)
- {
- if(ty<17)
- {
- z=cgetg(3,17);z[1]=lcopy(x);
- z[2]=lcopy(y);
- }
- else
- {
- if(ty!=19)
- {
- z=cgetg(ly+1,ty);z[1]=lcopy(x);
- for(i=2;i<=ly;i++)
- z[i]=lcopy(y[i-1]);
- }
- else
- {
- if(lg(y[1])!=2) err(concater);
- z=cgetg(ly+1,ty);p1=cgetg(2,18);
- z[1]=(long)p1;p1[1]=lcopy(x);
- for(i=2;i<=ly;i++)
- z[i]=lcopy(y[i-1]);
- }
- }
- }
- else
- {
- switch(tx)
- {
- case 17:
- if(ty<17)
- {
- z=cgetg(lx+1,tx);z[lx]=lcopy(y);
- for(i=1;i<lx;i++)
- z[i]=lcopy(x[i]);
- }
- else
- {
- switch(ty)
- {
- case 17: z=cgetg(lx+ly-1,tx);
- for(i=1;i<lx;i++)
- z[i]=lcopy(x[i]);
- for(i=1;i<ly;i++)
- z[lx+i-1]=lcopy(y[i]);
- break;
- case 18:
- if(lx==1) z=concat(x[1],y);
- else
- {
- if(ly!=1) err(concater);
- z=concat(x,y[1]);
- }
- break;
- case 19: if(lx!=ly) err(concater);
- z=cgetg(ly,ty);
- for(i=1;i<ly;i++)
- z[i]=lconcat(x[i],y[i]);
- break;
- default:;
- }
- }
- break;
- case 18:
- if(ty<17)
- {
- z=cgetg(lx+1,tx);z[lx]=lcopy(y);
- for(i=1;i<lx;i++)
- z[i]=lcopy(x[i]);
- }
- else
- {
- switch(ty)
- {
- case 17:
- if(lx==1) z=concat(x[1],y);
- else
- {
- if(ly!=1) err(concater);
- z=concat(x,y[1]);
- }
- break;
- case 18: z=cgetg(lx+ly-1,tx);
- for(i=1;i<lx;i++)
- z[i]=lcopy(x[i]);
- for(i=1;i<ly;i++)
- z[lx+i-1]=lcopy(y[i]);
- break;
- case 19: if(lx!=lg(y[1])) err(concater);
- z=cgetg(ly+1,ty);
- z[1]=lcopy (x);
- for(i=2;i<=ly;i++)
- z[i]=lcopy(y[i-1]);
- break;
- default:;
- }
- }
- break;
- case 19: dx=lg(x[1]);
- if(ty<17)
- {
- if(dx!=1) err(concater);
- z=cgetg(lx+1,tx);
- for(i=1;i<lx;i++)
- z[i]=lcopy(x[i]);
- p1=cgetg(2,18);z[lx]=(long)p1;
- p1[1]=lcopy(y);
- }
- else
- {
- switch(ty)
- {
- case 17: if(lx!=ly) err(concater);
- z=cgetg(lx,tx);
- for(i=1;i<lx;i++)
- z[i]=lconcat(x[i],y[i]);
- break;
- case 18: if(dx!=ly) err(concater);
- z=cgetg(lx+1,tx);
- for(i=1;i<lx;i++)
- z[i]=lcopy(x[i]);
- z[lx]=lcopy(y);
- break;
- case 19: if(dx!=lg(y[1])) err(concater);
- z=cgetg(lx+ly-1,tx);
- for(i=1;i<lx;i++)
- z[i]=lcopy(x[i]);
- for(i=1;i<ly;i++)
- z[lx+i-1]=lcopy(y[i]);
- break;
- default:;
- }
- }
- break;
- default:;
- }
- }
- return z;
- }
-
- GEN vectextract(x,l)
-
- GEN x,l;
-
- /* extraction des composantes de x suivants les bits du masque l */
- /* a usage interne donc aucune verification n'est faite. voir extract */
-
- {
- GEN y;
- long i,tx=typ(x),lx=lg(x),av,tetpil,f;
-
- if(!signe(l)) return cgetg(1,tx);
- av=avma;i=1;
- while(!mpodd(l))
- {
- l=shifti(l,-1);i++;
- }
- if(i>=lx) err(extracter3);
- l=shifti(l,-1);tetpil=avma;
- y=cgetg(2,tx);y[1]=lcopy(x[i]);
- i++;
- while(!gcmp0(l)&&(i<lx))
- {
- f=mpodd(l);l=shifti(l,-1);tetpil=avma;
- if(f) y=concat(y,x[i]);
- i++;
- }
- if(!gcmp0(l)) err(extracter3);
- y=gerepile(av,tetpil,y);
- return y;
- }
-
- GEN extract(x,l)
-
- GEN x,l;
-
- {
- long tl=typ(l),ll,lx,i,tx=typ(x),in;
- GEN y;
-
- if(tx<17) err(extracter1);
- if(tl==1) return vectextract(x,l);
- if((tl==17)||(tl==18))
- {
- ll=lg(l);y=cgetg(ll,tx);lx=lg(x);
- for(i=1;i<ll;i++)
- {
- in=itos(l[i]);if((in>=lx)||(in<=0)) err(extracter3);
- y[i]=lcopy(x[in]);
- }
- return y;
- }
- err(extracter2);
- }
-
- GEN matextract(x,l1,l2)
-
- GEN x,l1,l2;
-
- {
- GEN y;
- long av,tetpil;
-
- if(typ(x)!=19) err(matextrer);
- av=avma;y=extract(gtrans(extract(x,l2)),l1);tetpil=avma;
- return gerepile(av,tetpil,gtrans(y));
- }
-
-
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /* */
- /* OPERATIONS SCALAIRES-MATRICES */
- /* */
- /* ET DIVERS */
- /* */
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
-
-
- GEN gscalmat(x,n) /* cree la matrice carree n X n */
- /* contenant x*I */
- GEN x;
- long n;
-
- {
- long i,j,z;
- GEN y;
-
- z = lcopy(x);
- y=cgetg(n+1,19);
- for(i=1;i<=n;i++)
- {
- y[i]=lgetg(n+1,18);
- for(j=1;j<=n;j++)
- coeff(y,j,i)=(i==j ? z : zero);
- }
- return y;
- }
-
- GEN gscalsmat(x,n) /* idem au precedent avec x long du C */
-
- long x,n;
-
- {
- long i,j,z;
- GEN y;
-
- z=lstoi(x);
- y=cgetg(n+1,19);
- for(i=1;i<=n;i++)
- {
- y[i]=lgetg(n+1,18);
- for(j=1;j<=n;j++)
- coeff(y,j,i)=(i==j ? z : zero);
- }
- return y;
- }
-
- GEN idmat(n)
-
- long n;
- {
- return gscalmat(gun,n);
- }
-
- GEN gtomat(x)
-
- GEN x;
-
- {
- GEN y,p1;
- long tx=typ(x),lx,i;
-
- if(tx<17)
- {
- y=cgetg(2,19);p1=cgetg(2,18);y[1]=(long)p1;
- p1[1]=lcopy(x);
- }
- else
- switch(tx)
- {
- case 17: lx=lg(x);y=cgetg(lx,19);
- for(i=1;i<lx;i++)
- {
- p1=cgetg(2,18);p1[1]=lcopy(x[i]);
- y[i]=(long)p1;
- }
- break;
- case 18: y=cgetg(2,19);y[1]=lcopy(x);break;
- case 19: y=gcopy(x);break;
- }
- return y;
- }
-
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /* */
- /* ADDITION SCALAIRE + MATRICE */
- /* */
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
-
- GEN gaddmat(x,y) /* cree la matrice carree contenant x*I+y */
-
- GEN x,y;
-
- {
- long ly,dy,i,j;
- GEN z;
-
- ly=lg(y);dy=lg(y[1]);
- if((typ(y)!=19) || (ly!=dy)) err(gadmaer);
- z=cgetg(ly,19);
- for(i=1;i<ly;i++)
- {
- z[i]=lgetg(dy,18);
- for(j=1;j<dy;j++)
- coeff(z,j,i)=(i==j ? ladd(x,coeff(y,j,i)) : lcopy(coeff(y,j,i)));
- }
- return z;
- }
-
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /* */
- /* ADDITION SHORT + MATRICE */
- /* */
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
-
- GEN gaddsmat(s,y) /* idem au precedent avec x long du C */
-
- long s;
- GEN y;
-
- {
- long ly,dy,i,j;
- GEN z;
-
- ly=lg(y);dy=lg(y[1]);
- if((typ(y)!=19) || (ly!=dy)) err(gadsmaer);
- z=cgetg(ly,19);
- for(i=1;i<ly;i++)
- {
- z[i]=lgetg(dy,18);
- for(j=1;j<dy;j++)
- coeff(z,j,i)=(i==j ? laddsg(s,coeff(y,j,i)) : lcopy(coeff(y,j,i)));
- }
- return z;
- }
-
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /* */
- /* RESOLUTION DE A X=B */
- /* (METHODE DE GAUSS) */
- /* */
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
-
- GEN gauss(a,b)
-
- GEN a,b;
-
- {
- long p,u,m,nbli,nbco,i,j,k,av1,av2,av3,av4;
- GEN aa,x;
-
- if(typ(b)==19) return invmulmat(a,b);
- nbco=lg(a)-1;nbli=lg(a[1])-1;
- if (nbco!=nbli) err(gausser1);
- x=cgetg(nbli+1,18);av1=avma;
- for (j=1;j<=nbco;j++) x[j]=b[j];
- aa=cgetg(nbco+1,19);
- for (j=1;j<=nbco;j++)
- {
- aa[j]=lgetg(nbli+1,18);
- for (i=1;i<=nbli;i++) coeff(aa,i,j)=coeff(a,i,j);
- }
- for (i=1;i<nbli;i++)
- {
- p=coeff(aa,i,i);k=i;
- if (gcmp0(p))
- {
- for (k=i+1;(k<=nbli)&&gcmp0(coeff(aa,k,i));k++);
- if (k>nbco) err(matinv1);
- else
- {
- for (j=i;j<=nbco;j++)
- {
- u=coeff(aa,i,j);coeff(aa,i,j)=coeff(aa,k,j);
- coeff(aa,k,j)=u;
- }
- u=x[i];x[i]=x[k];x[k]=u;
- p=coeff(aa,i,i);
- }
- }
- for (k=i+1;k<=nbli;k++)
- {
- m=coeff(aa,k,i);
- if (!gcmp0(m))
- {
- m=ldiv(m,p);
- for (j=i+1;j<=nbco;j++)
- coeff(aa,k,j)=lsub(coeff(aa,k,j),gmul(m ,coeff(aa,i,j)));
- x[k]=lsub(x[k],gmul(m,x[i]));
- }
- }
- }
-
- /* Resolution systeme triangularise */
- av2=avma;
- p=coeff(aa,nbli,nbco);
- if (gcmp0(p)) err(matinv1);
- else
- {
- x[nbli]=ldiv(x[nbli],p);
- for (i=nbli-1;i>0;i--)
- {
- av3=avma;m=x[i];
- for (j=i+1;j<=nbco;j++)
- m= lsub(m,gmul(coeff(aa,i,j),x[j]));
- av4=avma;
- x[i]=lpile(av3,av4,gdiv(m,coeff(aa,i,i)));
- }
- }
- gerepile(av1,av2,1);
- return x;
- }
-
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /* */
- /* RANG D'UNE MATRICE m lignes x n colonnes */
- /* */
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
-
- long rank(x)
- GEN x;
-
- {
- GEN c,mun,p;
- long i,j,k,r,t,n,m,av;
-
- if (typ(x)!=19) err(kerer1);
- r=n=lg(x)-1;if(!r) return r;
- m=lg(x[1])-1;av=avma;
- x=gcopy(x);c=cgeti(m+1);
- mun=gneg(un);
- for(k=1;k<=m;k++) c[k]=0;
- for(k=1;k<=n;k++)
- {
- j=1;
- while((j<=m)&&(c[j]||gcmp0(coeff(x,j,k)))) j++;
- if (j<=m)
- {
- p=gdivsg(-1,coeff(x,j,k));
- coeff(x,j,k)=(long)mun;
- for(i=k+1;i<=n;i++) coeff(x,j,i)=lmul(p,coeff(x,j,i));
- for(t=1;t<=m;t++)
- if(t!=j)
- {
- p=(GEN)coeff(x,t,k);
- for(i=k+1;i<=n;i++) coeff(x,t,i)=ladd(coeff(x,t,i),gmul(p,coeff(x,j,i)));
- coeff(x,t,k)=zero;
- }
- c[j]=k;
- }
- else r--;
- }
- avma=av;
- return r;
- }
-
- GEN indexrank(x)
- GEN x;
-
- {
- GEN c,d,mun,p,y,p1,p2;
- long i,j,k,r,t,n,m,av,tetpil;
-
- if (typ(x)!=19) err(kerer1);
- r=n=lg(x)-1;if(!r) {y=cgetg(3,17);y[1]=lgetg(1,17);y[2]=lgetg(1,17);return y;}
- m=lg(x[1])-1;av=avma;
- x=gcopy(x);c=cgeti(m+1);d=cgeti(n+1);
- mun=gneg(un);
- for(k=1;k<=m;k++) c[k]=0;
- for(k=1;k<=n;k++)
- {
- j=1;
- while((j<=m)&&(c[j]||gcmp0(coeff(x,j,k)))) j++;
- if (j<=m)
- {
- p=gdivsg(-1,coeff(x,j,k));
- coeff(x,j,k)=(long)mun;
- for(i=k+1;i<=n;i++) coeff(x,j,i)=lmul(p,coeff(x,j,i));
- for(t=1;t<=m;t++)
- if(t!=j)
- {
- p=(GEN)coeff(x,t,k);
- for(i=k+1;i<=n;i++) coeff(x,t,i)=ladd(coeff(x,t,i),gmul(p,coeff(x,j,i)));
- coeff(x,t,k)=zero;
- }
- c[j]=k;d[k]=j;
- }
- else {r--;d[k]=0;}
- }
- p1=cgetg(r+1,17);p2=cgetg(r+1,17);
- for(i=0,k=1;k<=n;k++) if(d[k]) {p1[++i]=lstoi(d[k]);p2[i]=lstoi(k);}
- tetpil=avma;y=cgetg(3,17);y[1]=(long)sort(p1);y[2]=lcopy(p2);
- return gerepile(av,tetpil,y);
- }
-
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /* */
- /* NOYAU D'UNE MATRICE m lignes x n colonnes */
- /* ( Retourne une matrice de n-rang vecteurs independants ) */
- /* */
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
-
- GEN keri(x) /* Programme pour types ENTIERS */
- GEN x;
-
- {
- GEN c,d,y,v,pp;
- long i,j,k,r,t,n,n1,m,nbmot,av,av1;
- long p,p0,p1,q;
-
- if (typ(x)!=19) err(kerer1);
- n1=lg(x);n=n1-1;if(!n) return cgetg(1,19);
- m=lg(x[1])-1;av=avma;
- nbmot=200;
- c=cgetg(n1,19);
- for(j=1;j<=n;j++)
- {
- p=c[j]=lgetg(m+1,18);
- for(i=1;i<=m;i++)
- affii(coeff(x,i,j),((GEN)p)[i]=lgeti(nbmot));
- }
- x=c;
- p=un;
- pp=cgetg(n+1,18);
- for(j=1;j<=n;j++) pp[j]=lgeti(nbmot);
- c=cgeti(m+1);for(k=1;k<=m;k++) c[k]=0;
- d=cgeti(n1);
- av1=avma;
- for(r=0,k=1;k<=n;k++)
- {
- j=1;
- while((j<=m)&&(c[j]||!signe(coeff(x,j,k)))) j++;
- if (j<=m)
- {
- p0=p;p=coeff(x,j,k);
- for(t=1;t<=m;t++)
- if(t!=j)
- {
- q=coeff(x,t,k);
- for(i=k+1;i<=n;i++)
- {
- p1=lsubii(mulii(p,coeff(x,t,i)),mulii(q,coeff(x,j,i)));
- if(k>1) diviiz(p1,p0,coeff(x,t,i));
- else affii(p1,coeff(x,t,i));
- }
- }
- c[j]=k;d[k]=j;
- avma=av1;
- }
- else {r++;d[k]=0;affii(p,pp[k]);}
-
- }
- if(r) /* Il y a un noyau non nul */
- {
- av1=avma;
- y=cgetg(r+1,19);
- for(j=k=1;j<=r;j++,k++)
- {
- while(d[k]) k++;
- y[j]=(long)(v=cgetg(n1,18));
- for(i=1;i<k;i++) v[i]= d[i]? lcopy(coeff(x,d[i],k)) : zero;
- v[k]=lnegi(pp[k]);
- for(i=k+1;i<=n;i++) v[i]=zero;
- }
- return gerepile(av,av1,y);
- }
- else {avma=av;y=cgetg(1,19);return y;}
- }
-
- GEN ker(x) /* Programme pour types exacts */
- GEN x;
-
- {
- GEN c,d,y,mun,p;
- long i,j,k,r,t,n,n1,m,av,av1,av2;
-
- if (typ(x)!=19) err(kerer1);
- n1=lg(x);n=n1-1;if(!n) return cgetg(1,19);
- m=lg(x[1])-1;av=avma;x=gcopy(x);mun=gneg(un);r=0;
- c=cgeti(m+1);for(k=1;k<=m;k++) c[k]=0;
- d=cgeti(n1);
- av1=avma;
- for(k=1;k<=n;k++)
- {
- j=1;
- while((j<=m)&&(c[j]||gcmp0(coeff(x,j,k)))) j++;
- if (j<=m)
- {
-
- p=gdivsg(-1,coeff(x,j,k));
- coeff(x,j,k)=(long)mun;
- for(i=k+1;i<=n;i++) coeff(x,j,i)=lmul(p,coeff(x,j,i));
- for(t=1;t<=m;t++)
- if(t!=j)
- {
- p=(GEN)coeff(x,t,k);
- for(i=k+1;i<=n;i++) coeff(x,t,i)=ladd(coeff(x,t,i),gmul(p,coeff(x,j,i)));
- coeff(x,t,k)=zero;
- }
- c[j]=k;d[k]=j;
- av2=avma;
- x=gerepile(av1,av2,gcopy(x));
- }
- else {r++;d[k]=0;}
- }
- if(r)
- {
- av1=avma;
- y=cgetg(r+1,19);
- for(j=k=1;j<=r;j++,k++)
- {
- while(d[k]) k++;
- y[j]=(long)(p=cgetg(n1,18));
- for(i=1;i<k;i++) p[i]=d[i]? lcopy(coeff(x,d[i],k)):zero;
- p[k]=un;
- for(i=k+1;i<=n;i++) p[i]=zero;
- }
- return gerepile(av,av1,y);
- }
- else {avma=av;y=cgetg(1,19);return y;}
- }
-
- GEN image(x) /* Programme pour types exacts */
- GEN x;
-
- {
- GEN c,d,y,mun,p,x1;
- long i,j,k,r,t,n,n1,m,av,av1,av2;
-
- if (typ(x)!=19) err(kerer1);
- n1=lg(x);n=n1-1;if(!n) return cgetg(1,19);
- m=lg(x[1])-1;av=avma;x1=gcopy(x);mun=gneg(un);r=0;
- c=cgeti(m+1);for(k=1;k<=m;k++) c[k]=0;
- d=cgeti(n1);
- av1=avma;
- for(k=1;k<=n;k++)
- {
- j=1;
- while((j<=m)&&(c[j]||gcmp0(coeff(x1,j,k)))) j++;
- if (j<=m)
- {
-
- p=gdivsg(-1,coeff(x1,j,k));
- coeff(x1,j,k)=(long)mun;
- for(i=k+1;i<=n;i++) coeff(x1,j,i)=lmul(p,coeff(x1,j,i));
- for(t=1;t<=m;t++)
- if(t!=j)
- {
- p=(GEN)coeff(x1,t,k);
- for(i=k+1;i<=n;i++) coeff(x1,t,i)=ladd(coeff(x1,t,i),gmul(p,coeff(x1,j,i)));
- coeff(x1,t,k)=zero;
- }
- c[j]=k;d[k]=j;
- av2=avma;
- x1=gerepile(av1,av2,gcopy(x1));
- }
- else {r++;d[k]=0;}
- }
- if(r)
- {
- av1=avma;
- y=cgetg(n-r+1,19);
- for(j=k=1;j<=n-r;j++,k++)
- {
- while(!d[k]) k++;
- y[j]=lcopy(x[k]);
- }
- return gerepile(av,av1,y);
- }
- else {avma=av;return gcopy(x);}
- }
-
- GEN inverseimage(mat,y)
- GEN mat,y;
- {
- long av=avma,nbcol,i,j,l,tetpil;
- GEN met,noyau,lastcoeff,invimag;
-
- if ((typ(mat)!=19)||(typ(y)!=18)) err(kerer1);
- nbcol=lg(mat);
- met=cgetg(nbcol+1,19);
- for(j=1;j<=nbcol-1;j++) met[j]=mat[j];
- met[nbcol]=(long)y;
- noyau=ker(met);l=lg(noyau)-1;
- if(!l) {avma=av;return cgetg(1,18);}
- lastcoeff=gneg(coeff(noyau,nbcol,l));
- if(gcmp0(lastcoeff)) {avma=av;return cgetg(1,18);}
- tetpil=avma;
- invimag=cgetg(nbcol,18);
- for(i=1;i<=nbcol-1;i++) invimag[i]=ldiv(coeff(noyau,i,l),lastcoeff);
- return gerepile(av,tetpil,invimag);
- }
-
- GEN kerreel(x,prec) /* Programme pour types non exacts */
- /* gestion de pile a la fin seulement */
- GEN x;
- long prec;
- {
- GEN c,d,y,mun,p,eps;
- long i,j,k,r,t,n,n1,m,av,av1;
-
- if (typ(x)!=19) err(kerer1);
- n1=lg(x);n=n1-1;if(!n) return cgetg(1,19);
- m=lg(x[1])-1;av=avma;
- eps=cgetr(3);eps[2]=0x80000000;eps[1]=0x01800010-((prec-2)<<5);
- x=gcopy(x);
- mun=gneg(un);r=0;
- c=cgeti(m+1);for(k=1;k<=m;k++) c[k]=0;
- d=cgeti(n1);
- for(k=1;k<=n;k++)
- {
- j=1;
- while((j<=m)&&(c[j]||(gcmp(gabs(coeff(x,j,k),5),eps)<0))) j++;
- if (j<=m)
- {
-
- p=gdivsg(-1,coeff(x,j,k));
- coeff(x,j,k)=(long)mun;
- for(i=k+1;i<=n;i++) coeff(x,j,i)=lmul(p,coeff(x,j,i));
- for(t=1;t<=m;t++)
- if(t!=j)
- {
- p=(GEN)coeff(x,t,k);
- for(i=k+1;i<=n;i++) coeff(x,t,i)=ladd(coeff(x,t,i),gmul(p,coeff(x,j,i)));
- coeff(x,t,k)=zero;
- }
- c[j]=k;d[k]=j;
- }
- else{r++;d[k]=0;}
- }
- if(r)
- {
- av1=avma;
- y=cgetg(r+1,19);
- for(j=k=1;j<=r;j++,k++)
- {
- while(d[k]) k++;
- y[j]=(long)(p=cgetg(n1,18));
- for(i=1;i<k;i++) p[i]=d[i]? lcopy(coeff(x,d[i],k)):zero;
- p[k]=un;
- for(i=k+1;i<=n;i++) p[i]=zero;
- }
- return gerepile(av,av1,y);
- }
- else {avma=av;y=cgetg(1,19);return y;}
- }
-
- /* Etant donnee une matrice nxk de rang k<=n, on trouve une matrice nxn
- inversible dont les k premieres colonnes forment la matrice initiale;
- on ne verifie pas que les k colonnes sont lineairement independantes. */
-
- GEN suppl(x)
- GEN x;
- {
- long av=avma,tetpil,k,n,s,t;
- GEN y,p1,p2;
-
- if(typ(x)!=19) err(kerer1);
- k=lg(x)-1;if(!k) err(suppler1);
- n=lg(x[1])-1;if(k>n) err(suppler2);
- s=0;y=idmat(n);
- while(s<k)
- {
- s++;p1=gauss(y,x[s]);t=s;
- while((t<=n)&&gcmp0(p1[t])) t++;
- if(t>n) err(suppler2);
- p2=(GEN)y[s];y[s]=x[s];if(s!=t) y[t]=(long)p2;
- }
- tetpil=avma;return gerepile(av,tetpil,gcopy(y));
- }
-
- GEN image2(x)
- GEN x;
- {
- long av=avma,tetpil,k,n,i;
- GEN p1,p2;
-
- if(typ(x)!=19) err(kerer1);
- k=lg(x)-1;if(!k) return gcopy(x);
- n=lg(x[1])-1;p1=ker(x);k=lg(p1)-1;
- if(k) p1=suppl(p1);else p1=idmat(n);
- n=lg(p1)-1;
- tetpil=avma;p2=cgetg(n-k+1,19);
- for(i=k+1;i<=n;i++) p2[i-k]=lmul(x,p1[i]);
- return gerepile(av,tetpil,p2);
- }
-
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /* */
- /* VECTEURS PROPRES */
- /* (matrice de vecteurs propres independants */
- /* classes par valeurs propres croissantes ) */
- /* */
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
-
- GEN eigen(x,prec)
- GEN x;
- long prec;
-
- {
- GEN y,z,rr,p,ssesp,eps;
- long j,k,n,ly,av,av1,nbrac,nk,r1,r2,r3,flag;
-
-
- n=lg(x);
- av=avma;
- eps=cgetr(3);eps[2]=0x80000000;eps[1]=0x01800010-((prec-2)<<5);
- y=cgetg(n,19);ly=1;
- z=gcopy(x);
- p=caradj(x,0,0);rr=roots(p,prec);nbrac=lg(rr)-1;
- /* Bien sur ce n'est pas comme cela qu'on doit calculer les valeurs propres !*/
- for(k=1;k<=nbrac;k++)
- {
- r2=rr[k];flag=0;
- if(k>1) if(gcmp(gabs(gsub(r1,r2),5),eps)>0) flag=1;
- if(flag||(k==1))
- {
- r3=lround(r2);if(gcmp(gabs(gsub(r2,r3),5),eps)<0) r2=r3;
- {
- for(j=1;j<n;j++) coeff(z,j,j)=lsub(coeff(x,j,j),r2);
- ssesp=kerreel(z,prec);
- nk=lg(ssesp)-1;
- for(j=1;j<=nk;j++,ly++) y[ly]=ssesp[j];
- }
- r1=r2;
- }
- }
- z=cgetg(ly,19);
- av1=avma;
- for(k=1;k<ly;k++) z[k]=y[k];
- return gerepile(av,av1,gcopy(z));
- }
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /* */
- /* DETERMINANT */
- /* */
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
-
- /* ===================================================================*/
- /* Determinant types exacts : 1er pivot non nul */
- /*--------------------------------------------------------------------*/
-
- GEN det2(a)
-
- GEN a;
-
- {
- long p,u,m,nbli,nbco,i,j,k,av,av1,s;
- GEN aa,p1,x;
-
- if (typ(a)!=19) err(mattype1);
- nbco=lg(a)-1;nbli=lg(a[1])-1;
- if (nbco!=nbli) err(mattype1);
- av=avma;x=gun;s=1;
- aa=cgetg(nbco+1,19);
-
- for (j=1;j<=nbco;j++)
- {
- aa[j]=lgetg(nbli+1,18);
- for (i=1;i<=nbli;i++) coeff(aa,i,j)=coeff(a,i,j);
- }
-
- for (i=1;i<nbco;i++)
- {
- p=coeff(aa,i,i);k=i;
- if(gcmp0(p))
- {
-
- for (k=i+1;(k<=nbco)&&gcmp0(coeff(aa,i,k));k++);
- if (k>nbco)
- {
- avma=av;return gzero;
- }
- else
- {
- p=coeff(aa,i,k);
- u=aa[k];aa[k]=aa[i];aa[i]=u;
- s= -s;
- }
- }
- x=gmul(x,p);
-
- for (k=i+1;k<=nbco;k++)
- {
- m=coeff(aa,i,k);
- if (!gcmp0(m))
- {
- m=ldiv(m,p);
- for (j=i+1;j<=nbli;j++)
- {
- p1=gmul(m,coeff(aa,j,i));
- coeff(aa,j,k)=lsub(coeff(aa,j,k),p1);
- }
- }
- }
- }
- if(s<0) x=gneg(x);
- av1=avma;
- return gerepile(av,av1,gmul(x,coeff(aa,nbli,nbco)));
- }
-
- /* ===================================================================*/
- /* Determinant dans un anneau A : Tous les calculs dans A */
- /* division par le pivot precedent ( methode de Bareiss) */
- /*--------------------------------------------------------------------*/
-
- GEN det(a)
-
- GEN a;
-
- {
- long p,pprec,u,m,nbli,nbco,i,j,k,av,av1,s;
- GEN aa,p1;
-
-
- if (typ(a)!=19) err(mattype1);
- nbco=lg(a)-1;nbli=lg(a[1])-1;
- if (nbco!=nbli) err(mattype1);
- av=avma;
- aa=cgetg(nbco+1,19);
-
- for (j=1;j<=nbco;j++)
- {
- aa[j]=lgetg(nbli+1,18);
- for (i=1;i<=nbli;i++) coeff(aa,i,j)=coeff(a,i,j);
- }
- pprec=un;s=1;
- for (i=1;i<nbco;i++)
- {
- p=coeff(aa,i,i);k=i;
- if(gcmp0(p))
- {
-
- for (k=i+1;(k<=nbco)&&gcmp0(coeff(aa,i,k));k++);
- if (k>nbco)
- {
- avma=av;return gzero;
- }
- else
- {
- p=coeff(aa,i,k);
- u=aa[k];aa[k]=aa[i];aa[i]=u;
- s= -s;
- }
- }
- for (k=i+1;k<=nbco;k++)
- {
- m=coeff(aa,i,k);
- for (j=i+1;j<=nbli;j++)
- {
- p1=gsub(gmul(p,coeff(aa,j,k)),gmul(m,coeff(aa,j,i)));
- if((typ(p1)==10)&&(typ(pprec)==10)&&(varn(p1)==varn(pprec)))
- coeff(aa,j,k)=ldeuc(p1,pprec);
- else coeff(aa,j,k)=ldiv(p1,pprec);
- }
- }
- pprec=p;
- }
- av1=avma;
- return (s>0) ? gerepile(av,av1,gcopy(coeff(aa,nbli,nbco))) : gerepile(av,av1,gneg(coeff(aa,nbli,nbco)));
- }
-
- /* ===================================================================*/
- /* Determinant reel : pivot maximal */
- /*--------------------------------------------------------------------*/
-
-
- GEN detreel(a)
-
- GEN a;
-
- {
- long p,u,m,nbli,nbco,i,j,k,av,av1,s;
- GEN aa,p1,x;
-
- if (typ(a)!=19) err(mattype1);
- nbco=lg(a)-1;nbli=lg(a[1])-1;
- if (nbco!=nbli) err(mattype1);
- av=avma;s=1;x=gun;
- aa=cgetg(nbco+1,19);
-
- for (j=1;j<=nbco;j++)
- {
- aa[j]=lgetg(nbli+1,18);
- for (i=1;i<=nbli;i++) coeff(aa,i,j)=coeff(a,i,j);
- }
-
- for (i=1;i<nbco;i++)
- {
- p=labs(coeff(aa,i,i));k=i;
- for(j=i+1;j<=nbco;j++)
- if(gcmp(p1=gabs(coeff(aa,i,j),3),p)>0) {p=(long)p1;k=j;}
- if(gcmp0(p))
- {
- av1=av;return gerepile(av,av1,gcopy(p));
- }
- else
- {
- p=coeff(aa,i,k);
- if(k>i)
- {
- u=aa[k];aa[k]=aa[i];aa[i]=u;
- s= -s;
- }
- }
- x=gmul(x,p);
-
- for (k=i+1;k<=nbco;k++)
- {
- m=coeff(aa,i,k);
- if (!gcmp0(m))
- {
- m=ldiv(m,p);
- for (j=i+1;j<=nbli;j++)
- coeff(aa,j,k)=lsub(coeff(aa,j,k),gmul(m,coeff(aa,j,i)));
- }
- }
- }
- if(s<0) x=gneg(x);
- av1=avma;return gerepile(av,av1,gmul(x,coeff(aa,nbli,nbco)));
- }
-