home *** CD-ROM | disk | FTP | other *** search
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /* */
- /* BASE D'ENTIERS */
- /* */
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
-
- # include "genpari.h"
- GEN rquot(),ordmax(),rtran(),mtran(),matinv();
- void rowred();
-
- GEN allbase(x,code,y)
-
- GEN x,*y;
- long code;
-
- /*******************************************************************
- Entree: x polynome unitaire a coefficients dans Z de deg n
- definissant un corps de nombres K=Q(theta);
- code 0, 1 ou (long)p selon que l'on veut base, smallbase
- ou factoredbase.
- y pointeur sur un GEN destine a recevoir
- le discriminant du corps K.
- Sortie: retourne 1) un vecteur (horizontal) a n composantes,
- de polynomes a coeff dans Q (de deg 0,1...n-1)
- constituant une base de l'anneau des entiers de K.
- 2) le discriminant de K (dans *y).
- Rem: le denominateur commun des coef. est dans da.
- *******************************************************************/
-
- {
- GEN p,a,at,bt,b,da,db,q;
- long av=avma,tetpil,n,h,j,je,k,r,s,t,pro,v;
-
- if(typ(x)!=10) err(allbaser1);
- n=lgef(x)-3;if(n<=0) err(allbaser1);
- v=varn(x);*y=discsr(x);
- switch(code)
- {
- case 0: p=auxdecomp(absi(*y),1);h=lg(p[1])-1;break; /* base */
- case 1: p=auxdecomp(absi(*y),0);h=lg(p[1])-1;break; /* smallbase */
- default: p=(GEN)code;
- if((typ(p)!=19)||(lg(p)!=3)) err(factoreder1); /* factoredbase */
- h=lg(p[1])-1;
- q=gun;for(je=1;je<=h;je++) q=gmul(q,gpui(coeff(p,je,1),coeff(p,je,2)));
- if(gcmp(absi(q),absi(*y))) err(factoreder2);
- }
- a=idmat(n);da=gun;
- for(je=1;je<=h;je++)
- {
- if(gcmpgs(coeff(p,je,2),1)>0)
- {
- b=ordmax(x,coeff(p,je,1),coeff(p,je,2),&db);
- a=gmul(db,a);b=gmul(da,b);
- da=mulii(db,da);db=da;
- at=gtrans(a);bt=gtrans(b);
- for(r=n-1;r>=0;r--)
- for(s=r;s>=0;s--)
- while(signe(coef1(bt,s,r)))
- {
- q=rquot(coef1(at,s,s),coef1(bt,s,r));
- at[s+1]=(long)rtran(at[s+1],bt[r+1],q);
- for(t=s-1;t>=0;t--)
- {
- q=rquot(coef1(at,t,s),coef1(at,t,t));
- at[s+1]=(long)rtran(at[s+1],at[t+1],q);
- }
- pro=at[s+1];at[s+1]=bt[r+1];bt[r+1]=pro;
- }
- for(j=n-1;j>=0;j--)
- {
- for(k=0;k<j;k++)
- {
- while(signe(coef1(at,j,k)))
- {
- q=rquot(coef1(at,j,j),coef1(at,j,k));
- at[j+1]=(long)rtran(at[j+1],at[k+1],q);
- pro=at[j+1];at[j+1]=at[k+1];at[k+1]=pro;
- }
- }
- if(signe(coef1(at,j,j))<0)
- for(k=0;k<=j;k++) coef1(at,k,j)=lnegi(coef1(at,k,j));
- for(k=j+1;k<n;k++)
- {
- q=rquot(coef1(at,j,k),coef1(at,j,j));
- at[k+1]=(long)rtran(at[k+1],at[j+1],q);
- }
- }
- for(j=1;j<n;j++)
- if(!cmpii(coef1(at,j,j),coef1(at,j-1,j-1)))
- {
- coef1(at,0,j)=zero;
- for(k=1;k<=j;k++)
- coef1(at,k,j)=coef1(at,k-1,j-1);
- }
- a=gtrans(at);
- }
- }
- for(j=1;j<=n;j++)
- {
- *y=divii(mulii(coeff(a,j,j),*y),da);
- *y=divii(mulii(coeff(a,j,j),*y),da);
- }
- tetpil=avma;*y=gcopy(*y);at=cgetg(n+1,17);
- for(j=1;j<=n;j++)
- {
- q=cgetg(j+2,10);q[1]=0x1000002+j+(v<<16);at[j]=(long)q;
- for(k=2;k<=j+1;k++) q[k]=ldiv(coeff(a,j,k-1),da);
- }
- pro=lpile(av,tetpil,0)>>2;at+=pro;(*y)+=pro;
- return at;
- }
-
- GEN base(x,y)
- GEN x,*y;
- {
- return allbase(x,0,y);
- }
-
- GEN smallbase(x,y)
- GEN x,*y;
- {
- return allbase(x,1,y);
- }
-
- GEN factoredbase(x,p,y)
- GEN x,p,*y;
- {
- return allbase(x,(long)p,y);
- }
-
- GEN discf(x)
- GEN x;
- {
- GEN y;
- long av,tetpil;
-
- av=avma;allbase(x,0,&y);tetpil=avma;
- return gerepile(av,tetpil,gcopy(y));
- }
-
- GEN smalldiscf(x)
- GEN x;
- {
- GEN y;
- long av,tetpil;
-
- av=avma;allbase(x,1,&y);tetpil=avma;
- return gerepile(av,tetpil,gcopy(y));
- }
-
- GEN factoreddiscf(x,p)
- GEN x,p;
- {
- GEN y;
- long av,tetpil;
-
- av=avma;allbase(x,(long)p,&y);tetpil=avma;
- return gerepile(av,tetpil,gcopy(y));
- }
-
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /* */
- /* Quotient et Reste normalises ( -1/2 < r = x-q*y <= 1/2 ) */
- /* */
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
-
- GEN rquot(x,y)
-
- GEN x,y;
-
- {
- GEN u,v,w,p;
- long av,av1;
-
- av=avma;
- u=absi(y);v=shifti(x,1);w=shifti(y,1);
- if ( cmpii(u,v)>0) p=subii(v,u);
- else p=addsi(-1,addii(u,v));
- av1=avma;
- return(gerepile(av,av1,divii(p,w)));
- }
-
- GEN rrmdr(x,y)
-
- GEN x,y;
-
- {
- GEN p;
- long av,av1;
-
- av=avma;
- p=mulii(rquot(x,y),y);
- av1=avma;
- return(gerepile(av,av1,subii(x,p)));
- }
-
- GEN rinv(x,y)
-
- GEN x,y;
-
- {
- GEN a,c,q,r,t;
- long av,av1;
-
- av=avma;
- a=gun;c=gzero;
- while( signe(y))
- {
- q=rquot(x,y);
- r=subii(a,mulii(q,c));a=c;c=r;
- t=subii(x,mulii(q,y));x=y;y=t;
- }
- av1=avma;
- if (signe(x)<0) a=negi(a);
- if (signe(c)) { av1=avma; a=rrmdr(a,c); }
- return( gerepile(av,av1,a));
- }
-
- GEN rgcd(x,y)
-
- GEN x,y;
-
- {
- GEN z;
- long av,av1;
-
- av=avma;
- while(signe(y))
- {
- z=rrmdr(x,y);x=y;y=z;
- }
- av1=avma;
- return(gerepile(av,av1,absi(x)));
- }
-
-
- GEN rlcm(x,y)
-
- GEN x,y;
-
- {
- GEN d,z;
- long av,av1;
-
- av=avma;
- z=mulii(x,y);d=rgcd(x,y);
- av1=avma;
- return(gerepile(av,av1,divii(z,d)));
- }
-
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /* */
- /* Matrice compagnon du polynome unitaire x */
- /* */
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
-
- GEN companion(x)
-
- GEN x;
-
- {
- long i,j,l;
- GEN y;
-
- l=lgef(x)-2;y=cgetg(l,19);
- for(i=1;i<l;i++) y[i]=lgetg(l,18);
- for(i=0;i<l-2;i++)
- for(j=0;j<l-1;j++) coef1(y,i,j)=((i+1)==j) ? un : zero;
- for(j=0;j<l-1;j++) coef1(y,l-2,j)=lneg(x[j+2]);
- return(y);
- }
-
-
-
- GEN ordmax(f,p,e,ptdelta)
-
- GEN f,p,e;
- GEN *ptdelta;
-
- {
- GEN m,hh,pp,dd,ppdd,index,q,r,s,b,c,t,jp,v,delta;
- GEN cf[20],w[20],a;
- long h,i,j,k,sp,epsilon,n=lgef(f)-3,av=avma,tetpil,dec;
-
- a=cgetg(n*n+1,19);
- for(j=1;j<=n*n;j++)
- {
- a[j]=lgetg(n+1,18);
- for(k=1;k<=n;k++) coeff(a,k,j)=zero;
- }
- v=cgetg(n+1,18);
- cf[0]=idmat(n);
- cf[1]=companion(f);
- for(j=2;j<n;j++) cf[j]=gmul(cf[1],cf[j-1]);
- delta=gun; epsilon=itos(e);
- m=idmat(n);
-
- do
- {
- pp=mulii(p,p);
- dd=mulii(delta,delta);
- ppdd=mulii(dd,pp);
- b=matinv(m,delta);
- for(i=0;i<n;i++)
- {
- t=gscalsmat(0,n); /* t <--- matrice nulle d'ordre n */
- for(h=0;h<n;h++)
- for(j=0;j<n;j++)
- for(k=0;k<n;k++)
- coef1(t,j,k)=(long)rrmdr(addii(coef1(t,j,k),mulii(coef1(m,i,h),coef1(cf[h],j,k))),ppdd);
- c=gmul(t,b);
- w[i]=gmul(m,c);
- for(j=0;j<n;j++)
- for(k=0;k<n;k++)
- coef1(w[i],j,k)=(long)rrmdr(divii(coef1(w[i],j,k),dd),pp);
- }
- if(cmpis(p,n)>0)
- {
- for(i=0;i<n;i++)
- for(j=0;j<n;j++)
- {
- coeff(t,i+1,j+1)=zero;
- for(k=0;k<n;k++)
- for(h=0;h<n;h++)
- {
- r=modii(coef1(w[i],k,h),p);
- s=modii(coef1(w[j],h,k),p);
- coef1(t,i,j)=lmodii(addii(coef1(t,i,j),mulii(r,s)),p);
- }
- }
- }
- else
- {
- for(j=0;j<n;j++)
- {
- for(i=0;i<n;i++)
- coef1(b,i,j)=(i==0)? un:zero;
- /* ici la boucle en k calcule la puissance p mod p de w[j] */
- sp=itos(p);
- for(k=0;k<sp;k++)
- {
- for(i=0;i<n;i++)
- {
- v[i+1]=zero;
- for(h=0;h<n;h++)
- v[i+1]=lmodii(addii(v[i+1],mulii(coef1(b,h,j),coef1(w[j],h,i))),p);
- }
- for(i=0;i<n;i++) coef1(b,i,j)=v[i+1];
- }
- }
- q=p;t=b;
- while(cmpis(q,n)<0)
- {
- q=mulii(q,p);
- t=gmul(b,t);
- }
- }
- for(i=0;i<n;i++)
- for(j=0;j<n;j++)
- {
- coef1(a,j,i)=(i==j)? (long)p:zero;
- coef1(a,j,n+i)=lmodii(coef1(t,i,j),p);
- }
- rowred(a,2*n-1,pp);
- for(i=0;i<n;i++)
- for(j=0;j<n;j++)
- coef1(b,j,i)=coef1(a,j,i);
- jp=matinv(b,p);
- for(k=0;k<n;k++)
- {
- t=gmul(jp,w[k]);
- t=gmul(t,b);
- for(i=0;i<n;i++)
- for(j=0;j<n;j++)
- coef1(t,i,j)=ldivii(coef1(t,i,j),p);
- h=0;
- for(i=0;i<n;i++)
- for(j=0;j<n;j++)
- {
- coef1(a,k,h)=coef1(t,i,j);
- h++;
- }
- }
- rowred(a,n*n-1,pp);
- index=gun;
- for(i=0;i<n;i++)
- index=mulii(index,coef1(a,i,i));
- if (cmpsi(1,index))
- {
- delta=mulii(index,delta);
- for(i=0;i<n;i++)
- for(j=0;j<n;j++)
- coef1(c,i,j)=coef1(a,i,j);
- b=matinv(c,index);
- m=gmul(b,m);
- hh=delta;
- for(i=0;i<n;i++)
- for(j=0;j<n;j++)
- hh=rgcd(coef1(m,i,j),hh);
- if(cmpis(hh,1)>1)
- {
- delta=divii(delta,hh);
- for(i=0;i<n;i++)
- for(j=0;j<n;j++)
- coef1(m,i,j)=ldivii(coef1(m,i,j),hh);
- }
- q=index;
- while(!signe(modii(q,p)))
- {
- q=divii(q,p);
- epsilon=epsilon-2;
- }
- }
- }
- while(!gcmp1(index) && (epsilon>=2));
- tetpil=avma;delta=gcopy(delta);m=gcopy(m);
- dec=lpile(av,tetpil,0)>>2;
- *ptdelta=delta+dec;
- return m+dec;
- }
-
- void rowred(a,rlim,rmod)
- GEN a,rmod;
- long rlim;
-
- {
- long j,k,n,pro;
- GEN q;
-
- n=lg(a[1])-1;
- for(j=0;j<n;j++)
- {
- for(k=j+1;k<=rlim;k++)
- while (signe(coef1(a,j,k)))
- {
- q=rquot(coef1(a,j,j),coef1(a,j,k));
- a[j+1]=(long)mtran(a[j+1],a[k+1],q,rmod);
- pro=a[j+1];a[j+1]=a[k+1];a[k+1]=pro;
- }
- if (signe(coef1(a,j,j))<0)
- for(k=j;k<n;k++) coef1(a,k,j)=lnegi(coef1(a,k,j));
- for(k=0;k<j;k++)
- {
- q=rquot(coef1(a,j,k),coef1(a,j,j));
- a[k+1]=(long)mtran(a[k+1],a[j+1],q,rmod);
- }
- }
- }
-
- GEN rtran(v,w,q)
- GEN v,w,q ;
-
- {
- long av,tetpil;
- GEN p1;
-
- if (signe(q))
- {
- av=avma;p1=gmul(q,w);tetpil=avma;
- return gerepile(av,tetpil,gsub(v,p1));
- }
- else return v;
- }
-
- GEN mtran(v,w,q,m)
- GEN v,w,q,m;
-
- {
- long k;
-
- if (signe(q))
- {
- for(k=0;k<lg(v)-1;k++)
- {
- v[k+1]=(long)rrmdr(subii(v[k+1],modii(mulii(q,w[k+1]),m)),m);
- }
- }
- return v;
- }
-
-
- GEN matinv(x,d)
- GEN x,d;
- /*=======================================================================
- Calcule d/x ou d est entier et x matrice triangulaire inferieure
- entiere dont les coeff diagonaux divisent d ( resultat entier).
- ========================================================================*/
- {
- long n,h,i,j,k,av,av1;
- GEN y;
-
- av=avma;
- y=idmat(n=lg(x)-1);
- for(i=1;i<=n;i++)
- coeff(y,i,i)=ldivii(d,coeff(x,i,i));
- for(i=2;i<=n;i++)
- for(j=i-1;j;j--)
- {
- for(h=zero,k=j+1;k<=i;k++)
- h=ladd(h,mulii(coeff(y,i,k),coeff(x,k,j)));
- coeff(y,i,j)=ldivii(negi(h),coeff(x,j,j));
- }
- av1=avma;
- return gerepile(av,av1,gcopy(y));
- }
-
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~ ~*/
- /*~ HERMITE NORMAL FORM REDUCTION ~*/
- /*~ ~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
-
- GEN hnfold(x)
- GEN x;
- {
- long li,co,av,tetpil,i,j,j1,def,ind,lim;
- GEN p1,p2,y,z,m1;
-
- if(typ(x)!=19) err(hnfer1);
- lim=(avma+bot)>>1;
- av=avma;co=lg(x);li=lg(x[1]);y=gcopy(x);
- def=co;
- if(li>co) err(hnfer2);
- for(i=li-1;i>=1;i--)
- {
- def--;j=1;while((j<def)&&(!signe(coeff(y,i,j)))) j++;
- while(j<def)
- {
- m1=absi(coeff(y,i,j));ind=j;
- for(j1=j+1;j1<def;j1++)
- {
- p1=(GEN)coeff(y,i,j1);
- if(signe(p1)&&(cmpii(absi(p1),m1)<0))
- {
- m1=p1;ind=j1;
- }
- }
- p1=(GEN)coeff(y,i,def);
- if(!(signe(p1)&&(cmpii(absi(p1),m1)<0)))
- {
- p1=(GEN)y[def];y[def]=y[ind];y[ind]=(long)p1;
- }
- p1=(GEN)coeff(y,i,def);if(signe(p1)<0) y[def]=lneg(y[def]);
- p1=(GEN)coeff(y,i,def);
- for(j=1;j<def;j++)
- {
- p2=gdivent(coeff(y,i,j),p1);
- y[j]=lsub(y[j],gmul(p2,y[def]));
- }
- j=1;while((j<def)&&(!signe(coeff(y,i,j)))) j++;
- }
- p1=(GEN)coeff(y,i,def);
- if(signe(p1)<0) {y[def]=lneg(y[def]);p1=(GEN)coeff(y,i,def);}
- if(signe(p1))
- {
- for(j=def+1;j<co;j++)
- {
- p2=gdivent(coeff(y,i,j),p1);
- y[j]=lsub(y[j],gmul(p2,y[def]));
- }
- }
- if(avma<lim) {tetpil=avma;y=gerepile(av,tetpil,gcopy(y));}
- }
- tetpil=avma;z=cgetg(li,19);
- for(j=1;j<li;j++)
- {
- z[j]=lcopy(y[j+co-li]);
- }
- for(j=1;j<=co-li;j++) if(!gcmp0(y[j])) err(hnfer3);
- return gerepile(av,tetpil,z);
- }
-
- GEN hnfold2(x)
- GEN x;
- {
- long li,co,av,tetpil,i,j,def,lim;
- GEN p1,p2,y,z,u,v,d;
-
- if(typ(x)!=19) err(hnfer1);
- lim=(avma+bot)>>1;
- av=avma;co=lg(x);li=lg(x[1]);y=gcopy(x);
- def=co;
- if(li>co) err(hnfer2);
- for(i=li-1;i>=1;i--)
- {
- def--;j=def;while(j&&(!signe(coeff(y,i,j)))) j--;
- while(j>1)
- {
- d=bezout(coeff(y,i,j),coeff(y,i,j-1),&u,&v);
- p1=gadd(gmul(u,y[j]),gmul(v,y[j-1]));
- y[j]=lsub(gmul(divii(coeff(y,i,j),d),y[j-1]),gmul(divii(coeff(y,i,j-1),d),y[j]));
- y[j-1]=(long)p1;
- j--;while(j&&(!signe(coeff(y,i,j)))) j--;
- }
- if(j==1)
- {
- d=bezout(coeff(y,i,1),coeff(y,i,def),&u,&v);
- p1=gadd(gmul(u,y[1]),gmul(v,y[def]));
- y[1]=lsub(gmul(divii(coeff(y,i,1),d),y[def]),gmul(divii(coeff(y,i,def),d),y[1]));
- y[def]=(long)p1;
- }
- p1=(GEN)coeff(y,i,def);
- if(signe(p1)<0) {y[def]=lneg(y[def]);p1=(GEN)coeff(y,i,def);}
- if(signe(p1))
- {
- for(j=def+1;j<co;j++)
- {
- p2=gdivent(coeff(y,i,j),p1);
- y[j]=lsub(y[j],gmul(p2,y[def]));
- }
- }
- if(avma<lim) {tetpil=avma;y=gerepile(av,tetpil,gcopy(y));}
- }
- tetpil=avma;z=cgetg(li,19);
- for(j=1;j<li;j++)
- {
- z[j]=lcopy(y[j+co-li]);
- }
- for(j=1;j<=co-li;j++) if(!gcmp0(y[j])) err(hnfer3);
- return gerepile(av,tetpil,z);
- }
-
- GEN hnf(x)
- GEN x;
- {
- long li,co,av,tetpil,i,j,k,def,ldef,lim;
- GEN p1,p2,y,z,u,v,d;
-
- if(typ(x)!=19) err(hnfer1);
- lim=(avma+bot)>>1;
- av=avma;co=lg(x);li=lg(x[1]);y=gcopy(x);
- def=co;ldef=(li>co)?li-co+1:1;
- for(i=li-1;i>=ldef;i--)
- {
- def--;j=def-1;while(j&&(!signe(coeff(y,i,j)))) j--;
- while(j>1)
- {
- d=bezout(coeff(y,i,j),coeff(y,i,j-1),&u,&v);
- p1=gadd(gmul(u,y[j]),gmul(v,y[j-1]));
- y[j]=lsub(gmul(divii(coeff(y,i,j),d),y[j-1]),gmul(divii(coeff(y,i,j-1),d),y[j]));
- y[j-1]=(long)p1;
- j--;while(j&&(!signe(coeff(y,i,j)))) j--;
- }
- if(j==1)
- {
- d=bezout(coeff(y,i,1),coeff(y,i,def),&u,&v);
- p1=gadd(gmul(u,y[1]),gmul(v,y[def]));
- y[1]=lsub(gmul(divii(coeff(y,i,1),d),y[def]),gmul(divii(coeff(y,i,def),d),y[1]));
- y[def]=(long)p1;
- }
- p1=(GEN)coeff(y,i,def);
- if(signe(p1)<0) {y[def]=lneg(y[def]);p1=(GEN)coeff(y,i,def);}
- if(signe(p1))
- {
- for(j=def+1;j<co;j++)
- {
- p2=gdivent(coeff(y,i,j),p1);
- y[j]=lsub(y[j],gmul(p2,y[def]));
- }
- }
- else def++;
- if(avma<lim) {tetpil=avma;y=gerepile(av,tetpil,gcopy(y));}
- }
- for(i=0,j=1;j<co;j++) if(!gcmp0(y[j])) i++;
- tetpil=avma;z=cgetg(i+1,19);
- for(k=0,j=1;j<co;j++) if(!gcmp0(y[j])) z[++k]=lcopy(y[j]);
- return gerepile(av,tetpil,z);
- }
-
-
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~ ~*/
- /*~ SMITH NORMAL FORM REDUCTION ~*/
- /*~ ~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
- /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
-
- GEN smith2(x)
- GEN x;
- {
- long li,av,tetpil,i,j,k,l,lim,c,fl,n;
- GEN p1,p2,p3,p4,y,z,b,u,v,d;
-
- if(typ(x)!=19) err(hnfer1);
- lim=(avma+bot)>>1;
- av=avma;n=lg(x)-1;li=lg(x[1])-1;y=gcopy(x);
- if(li!=n) err(hnfer2);
- for(i=n;i>=2;i--)
- {
- do
- {
- for(j=i-1;j>=1;j--)
- {
- p1=(GEN)coeff(y,i,j);
- if(signe(p1))
- {
- p2=(GEN)coeff(y,i,i);
- if(gegal(p1,p2)) {d=p1;u=gun;v=gzero;p3=gun;p4=gun;}
- else if(!signe(addii(p1,p2)))
- {
- d=absi(p1);u=(signe(p2)>0)?gun:gneg(un);
- v=gzero;p3=u;p4=gneg(u);
- }
- else {d=bezout(p2,p1,&u,&v);p3=divii(p2,d);p4=divii(p1,d);}
- for(k=1;k<=i;k++)
- {
- b=addii(mulii(u,coeff(y,k,i)),mulii(v,coeff(y,k,j)));
- coeff(y,k,j)=lsubii(mulii(p3,coeff(y,k,j)),mulii(p4,coeff(y,k,i)));
- coeff(y,k,i)=(long)b;
- }
- }
- }
- c=0;
- for(j=i-1;j>=1;j--)
- {
- p1=(GEN)coeff(y,j,i);
- if(signe(p1))
- {
- p2=(GEN)coeff(y,i,i);
- if(gegal(p1,p2)) {d=p1;u=gun;v=gzero;p3=gun;p4=gun;}
- else if(!signe(addii(p1,p2)))
- {
- d=absi(p1);u=(signe(p2)>0)?gun:gneg(un);
- v=gzero;p3=u;p4=gneg(u);
- }
- else {d=bezout(p2,p1,&u,&v);p3=divii(p2,d);p4=divii(p1,d);}
- for(k=1;k<=i;k++)
- {
- b=addii(mulii(u,coeff(y,i,k)),mulii(v,coeff(y,j,k)));
- coeff(y,j,k)=lsubii(mulii(p3,coeff(y,j,k)),mulii(p4,coeff(y,i,k)));
- coeff(y,i,k)=(long)b;
- }
- c++;
- }
- }
- if(!c)
- {
- b=(GEN)coeff(y,i,i);fl=1;
- if(signe(b))
- {
- for(k=1;(k<i)&&fl;k++)
- for(l=1;(l<i)&&fl;l++)
- fl= !signe(modii(coeff(y,k,l),b));
- }
- if(!fl) {l--;y[i]=ladd(y[i],y[l]);}
- }
- if(avma<lim) {tetpil=avma;y=gerepile(av,tetpil,gcopy(y));}
- }
- while(c||(!fl));
- }
- tetpil=avma;z=cgetg(n+1,17);
- for(k=1;k<=n;k++) z[k]=labs(coeff(y,k,k));
- return gerepile(av,tetpil,z);
- }
-
- GEN smith(x)
- GEN x;
- {
- long li,av,tetpil,i,j,k,l,lim,c,fl,n;
- GEN p1,p2,p3,p4,y,z,b,u,v,d;
-
- if(typ(x)!=19) err(hnfer1);
- lim=(avma+bot)>>1;
- av=avma;n=lg(x)-1;if(!n) return cgetg(1,17);
- li=lg(x[1])-1;y=gcopy(x);
- if(li!=n) err(hnfer2);
- for(i=n;i>=2;i--)
- {
- do
- {
- c=0;
- for(j=i-1;j>=1;j--)
- {
- p1=(GEN)coeff(y,i,j);
- if(signe(p1))
- {
- p2=(GEN)coeff(y,i,i);
- if(gegal(p1,p2)) {d=p1;u=gun;v=gzero;p3=gun;p4=gun;}
- else if(!signe(addii(p1,p2)))
- {
- d=absi(p1);u=(signe(p2)>0)?gun:gneg(un);
- v=gzero;p3=u;p4=gneg(u);
- }
- else {d=bezout(p2,p1,&u,&v);p3=divii(p2,d);p4=divii(p1,d);}
- for(k=1;k<=i;k++)
- {
- b=addii(mulii(u,coeff(y,k,i)),mulii(v,coeff(y,k,j)));
- coeff(y,k,j)=lsubii(mulii(p3,coeff(y,k,j)),mulii(p4,coeff(y,k,i)));
- coeff(y,k,i)=(long)b;
- }
- c++;
- }
- }
- for(j=i-1;j>=1;j--)
- {
- p1=(GEN)coeff(y,j,i);
- if(signe(p1))
- {
- p2=(GEN)coeff(y,i,i);
- if(gegal(p1,p2)) {d=p1;u=gun;v=gzero;p3=gun;p4=gun;}
- else if(!signe(addii(p1,p2)))
- {
- d=absi(p1);u=(signe(p2)>0)?gun:gneg(un);
- v=gzero;p3=u;p4=gneg(u);
- }
- else {d=bezout(p2,p1,&u,&v);p3=divii(p2,d);p4=divii(p1,d);}
- for(k=1;k<=i;k++)
- {
- b=addii(mulii(u,coeff(y,i,k)),mulii(v,coeff(y,j,k)));
- coeff(y,j,k)=lsubii(mulii(p3,coeff(y,j,k)),mulii(p4,coeff(y,i,k)));
- coeff(y,i,k)=(long)b;
- }
- }
- }
- if(!c)
- {
- b=(GEN)coeff(y,i,i);fl=1;
- if(signe(b))
- {
- for(k=1;(k<i)&&fl;k++)
- for(l=1;(l<i)&&fl;l++)
- fl= !signe(modii(coeff(y,k,l),b));
- }
- if(!fl)
- {
- k--;
- for(l=1;l<=i;l++)
- coeff(y,i,l)=laddii(coeff(y,i,l),coeff(y,k,l));
- }
- }
- if(avma<lim) {tetpil=avma;y=gerepile(av,tetpil,gcopy(y));}
- }
- while(c||(!fl));
- }
- tetpil=avma;z=cgetg(n+1,17);
- for(j=0,k=1;k<=n;k++) if(!signe(coeff(y,k,k))) z[++j]=zero;
- for(k=1;k<=n;k++) if(signe(coeff(y,k,k))) z[++j]=labs(coeff(y,k,k));
- return gerepile(av,tetpil,z);
- }
-
-
- GEN oldsmith(x)
- GEN x;
- {
- long li,av,tetpil,i,j,k,l,lim,c,fl,n;
- GEN p1,p2,p3,p4,y,z,b,u,v,d;
-
- if(typ(x)!=19) err(hnfer1);
- lim=(avma+bot)>>1;
- av=avma;n=lg(x)-1;li=lg(x[1])-1;y=gcopy(x);
- if(li!=n) err(hnfer2);
- for(i=n;i>=2;i--)
- {
- do
- {
- c=0;
- for(j=i-1;j>=1;j--)
- {
- p1=(GEN)coeff(y,i,j);
- if(signe(p1))
- {
- p2=(GEN)coeff(y,i,i);
- if(gegal(p1,p2)) {d=p1;u=gun;v=gzero;p3=gun;p4=gun;}
- else {d=bezout(p2,p1,&u,&v);p3=divii(p2,d);p4=divii(p1,d);}
- for(k=1;k<=i;k++)
- {
- b=addii(mulii(u,coeff(y,k,i)),mulii(v,coeff(y,k,j)));
- coeff(y,k,j)=lsubii(mulii(p3,coeff(y,k,j)),mulii(p4,coeff(y,k,i)));
- coeff(y,k,i)=(long)b;
- }
- c++;
- }
- }
- for(j=i-1;j>=1;j--)
- {
- p1=(GEN)coeff(y,j,i);
- if(signe(p1))
- {
- p2=(GEN)coeff(y,i,i);
- if(gegal(p1,p2)) {d=p1;u=gun;v=gzero;p3=gun;p4=gun;}
- else {d=bezout(p2,p1,&u,&v);p3=divii(p2,d);p4=divii(p1,d);}
- for(k=1;k<=i;k++)
- {
- b=addii(mulii(u,coeff(y,i,k)),mulii(v,coeff(y,j,k)));
- coeff(y,j,k)=lsubii(mulii(p3,coeff(y,j,k)),mulii(p4,coeff(y,i,k)));
- coeff(y,i,k)=(long)b;
- }
- }
- }
- if(!c)
- {
- b=(GEN)coeff(y,i,i);fl=1;
- if(signe(b))
- {
- for(k=1;(k<i)&&fl;k++)
- for(l=1;(l<i)&&fl;l++)
- fl= !signe(modii(coeff(y,k,l),b));
- }
- if(!fl)
- {
- k--;
- for(l=1;l<=i;l++)
- coeff(y,i,l)=laddii(coeff(y,i,l),coeff(y,k,l));
- }
- }
- if(avma<lim) {tetpil=avma;y=gerepile(av,tetpil,gcopy(y));}
- }
- while(c||(!fl));
- }
- tetpil=avma;z=cgetg(n+1,17);
- for(j=0,k=1;k<=n;k++) if(!signe(coeff(y,k,k))) z[++j]=zero;
- for(k=1;k<=n;k++) if(signe(coeff(y,k,k))) z[++j]=labs(coeff(y,k,k));
- return gerepile(av,tetpil,z);
- }
-
-
- GEN transroot(x,i,j)
- GEN x;
- int i,j;
- {
- long n=lg(x),k;
- GEN y;
-
- y=cgetg(n,18);
- for(k=1;k<n;k++) y[k]=((k==i)||(k==j))?x[i+j-k]:x[k];
- return y;
- }
-
- GEN tchirnhausen(x)
- GEN x;
- {
- long av=avma,tetpil,v,n,a,b,c;
- GEN u;
-
- if(typ(x)!=10) err(galer1);
- n=lgef(x)-3;if(n<=0) err(galer1);v=varn(x);
- if(v) {u=gcopy(x);setvarn(u,0);x=u;}
- do
- {
- a=rand()&15;if(!a) a=1;b=rand()&31;if(b>=16) b-=32;
- c=rand()&31;if(c>=16) c-=32;
- u=caract(gmodulcp(gaddsg(c,gmul(polx[0],gaddsg(b,gmulsg(a,polx[0])))),x),v);
- }
- while(lgef(polgcd(u,deriv(u,v)))>=4);
- tetpil=avma;return gerepile(av,tetpil,gcopy(u));
- }
-
- int gpolcomp(p1,p2)
- GEN p1,p2;
- {
- int d,j;
-
- d=lgef(p1)-3;if((lgef(p2)-3)!=d) err(gpolcompbug1);
- j=d+1;while((j>=2)&&gegal(absi(p1[j]),absi(p2[j]))) j--;
- if(j==1) return 0;
- return gcmp(absi(p1[j]),absi(p2[j]));
- }
-
- GEN galois(x,prec)
- GEN x;
- long prec;
- {
-
- long av=avma,av1,i,j,k,n,f,l,l2,e,e1;
- GEN x1,p1,p2,p3,p4,p5,p6,y,z;
- static int ind5[20]={2,5,3,4,1,3,4,5,1,5,2,4,1,2,3,5,1,4,2,3};
- static int ind6[60]={3,5,4,6,2,6,4,5,2,3,5,6,2,4,3,6,2,5,3,4,1,4,5,6,1,5,3,6,1,6,3,4,1,3,4,5,1,6,2,5,1,2,4,6,1,5,2,4,1,3,2,6,1,2,3,5,1,4,2,3};
-
- if(typ(x)!=10) err(galer1);
- n=lgef(x)-3;if(n<=0) err(galer1);
- if(n>7) err(impl,"galois of degree higher than 7");
- x=gdiv(x,content(x));
- for(i=2;i<=n+2;i++) if(typ(x[i])!=1) err(galer2);
- p1=(GEN)x[n+2];
- if(!gcmp1(p1))
- {
- x1=cgetg(n+3,10);x1[1]=x[1];x1[n+2]=un;p2=gun;
- for(i=n+1;i>=2;i--) {x1[i]=lmul(x[i],p2);if(i>2) p2=gmul(p1,p2);}
- x=x1;
- }
- p1=factor(x);
- if((lg(p1[1])>2)||(!gcmp1(coeff(p1,1,2)))) err(impl,"galois of reducible polynomial");
- x1=gcopy(x);av1=avma;
- for(;;)
- {
- switch(n)
- {
- case 1: avma=av;y=cgetg(4,17);y[1]=y[3]=un;y[2]=lneg(gun);return y;
- case 2: avma=av;y=cgetg(4,17);y[1]=deux;y[3]=un;y[2]=lneg(gun);return y;
- case 3: f=carreparfait(discsr(x));avma=av;y=cgetg(4,17);y[3]=un;
- if(f) {y[2]=un;y[1]=lstoi(3);return y;}
- else {y[2]=lneg(gun);y[1]=lstoi(6);return y;}
- case 4:
- do
- {
- p1=roots(x,prec);p2=p1;
- p3=gzero;for(i=1;i<=4;i++) p3=gadd(p3,gmul(p2[i],gsqr(p2[(i&3)+1])));
- p4=gsub(polx[0],p3);p2=transroot(p1,1,2);
- p3=gzero;for(i=1;i<=4;i++) p3=gadd(p3,gmul(p2[i],gsqr(p2[(i&3)+1])));
- p4=gmul(p4,gsub(polx[0],p3));p2=transroot(p1,1,3);
- p3=gzero;for(i=1;i<=4;i++) p3=gadd(p3,gmul(p2[i],gsqr(p2[(i&3)+1])));
- p4=gmul(p4,gsub(polx[0],p3));p2=transroot(p1,1,4);
- p3=gzero;for(i=1;i<=4;i++) p3=gadd(p3,gmul(p2[i],gsqr(p2[(i&3)+1])));
- p4=gmul(p4,gsub(polx[0],p3));p2=transroot(p1,2,3);
- p3=gzero;for(i=1;i<=4;i++) p3=gadd(p3,gmul(p2[i],gsqr(p2[(i&3)+1])));
- p4=gmul(p4,gsub(polx[0],p3));p2=transroot(p1,3,4);
- p3=gzero;for(i=1;i<=4;i++) p3=gadd(p3,gmul(p2[i],gsqr(p2[(i&3)+1])));
- p4=gmul(p4,gsub(polx[0],p3));p5=grndtoi(greal(p4),&e);
- e=max(e,gexpo(gimag(p4)));
- if(e> -10) prec=(prec<<2)-2;
- }
- while(e> -10);
- p6=ggcd(p5,deriv(p5,0));
- f=(typ(p6)==10)&&(lgef(p6)>3);
- if(f) goto tchi;
- p1=factor(p5);p2=(GEN)p1[1];l=lg(p2)-1;
- switch(l)
- {
- case 1: f=carreparfait(discsr(x));avma=av;y=cgetg(4,17);y[3]=un;
- if(f) {y[2]=un;y[1]=lstoi(12);return y;}
- else {y[2]=lneg(gun);y[1]=lstoi(24);return y;}
- case 2: avma=av;y=cgetg(4,17);y[3]=un;y[2]=lneg(gun);y[1]=lstoi(8);return y;
- case 3: l2=lgef(p2[1])-3;
- if(l2==2) {avma=av;y=cgetg(4,17);y[3]=y[2]=un;y[1]=lstoi(4);return y;}
- else {avma=av;y=cgetg(4,17);y[3]=un;y[2]=lneg(gun);y[1]=lstoi(4);return y;}
- default: err(galbug1);
- }
- case 5:
- do
- {
- do
- {
- p1=roots(x,prec);z=cgetg(7,17);
- for(l=1;l<=5;l++)
- {
- p2=(l==1)?p1:transroot(p1,1,l);
- p3=gzero;k=0;for(i=1;i<=5;i++)
- {
- p5=gadd(gmul(p2[ind5[k]],p2[ind5[k+1]]),gmul(p2[ind5[k+2]],p2[ind5[k+3]]));
- p3=gadd(p3,gmul(gsqr(p2[i]),p5));k+=4;
- }
- z[l]=lrndtoi(greal(p3),&e);
- p4=(l==1)?gsub(polx[0],p3):gmul(p4,gsub(polx[0],p3));
- }
- p2=transroot(p1,2,5);
- p3=gzero;k=0;for(i=1;i<=5;i++)
- {
- p5=gadd(gmul(p2[ind5[k]],p2[ind5[k+1]]),gmul(p2[ind5[k+2]],p2[ind5[k+3]]));
- p3=gadd(p3,gmul(gsqr(p2[i]),p5));k+=4;
- }
- z[6]=lrndtoi(greal(p3),&e);
- p4=gmul(p4,gsub(polx[0],p3));
- p5=grndtoi(greal(p4),&e);
- e=max(e,gexpo(gimag(p4)));
- if(e> -10) prec=(prec<<2)-2;
- }
- while(e> -10);
- p6=ggcd(p5,deriv(p5,0));
- f=(typ(p6)==10)&&(lgef(p6)>3);
- if(f) goto tchi;
- p3=factor(p5);l=lg(p3[1])-1;
- f=carreparfait(discsr(x));
- if(l==1)
- {
- avma=av;y=cgetg(4,17);y[3]=un;
- if(f) {y[2]=un;y[1]=lstoi(60);return y;}
- else {y[2]=lneg(gun);y[1]=lstoi(120);return y;}
- }
- else
- {
- if(f)
- {
- l=1;while((l<=6)&&(!gcmp0(poleval(p5,z[l])))) l++;
- if(l>6) err(galbug4);
- p2=(l==6)?transroot(p1,2,5):transroot(p1,1,l);
- p3=gzero;
- for(i=1;i<=5;i++)
- {
- j=(i%5)+1;
- p3=gadd(p3,gmul(gmul(p2[i],p2[j]),gsub(p2[j],p2[i])));
- }
- p5=gmul(p3,p3);p4=grndtoi(greal(p5),&e1);
- e1=max(e1,gexpo(gimag(p5)));
- if(e1<= -10)
- {
- if(gcmp0(p4)) goto tchi;
- f=carreparfait(p4);
- avma=av;y=cgetg(4,17);y[3]=y[2]=un;y[1]=lstoi(f?5:10);return y;
- }
- else prec=(prec<<2)-2;
- }
- else
- {
- avma=av;y=cgetg(4,17);y[3]=un;y[2]=lneg(gun);y[1]=lstoi(20);return y;
- }
- }
- }
- while(e1> -10);
- case 6:
- do
- {
- do
- {
- p1=roots(x,prec);
- for(l=1;l<=6;l++)
- {
- p2=(l==1)?p1:transroot(p1,1,l);
- p3=gzero;k=0;for(i=1;i<=5;i++) for(j=i+1;j<=6;j++)
- {
- p5=gadd(gmul(p2[ind6[k]],p2[ind6[k+1]]),gmul(p2[ind6[k+2]],p2[ind6[k+3]]));
- p3=gadd(p3,gmul(gsqr(gmul(p2[i],p2[j])),p5));k+=4;
- }
- p4=(l==1)?gsub(polx[0],p3):gmul(p4,gsub(polx[0],p3));
- }
- p5=grndtoi(greal(p4),&e);
- e=max(e,gexpo(gimag(p4)));
- if(e> -10) prec=(prec<<2)-2;
- }
- while(e> -10);
- p6=ggcd(p5,deriv(p5,0));
- f=(typ(p6)==10)&&(lgef(p6)>3);
- if(f) goto tchi;
- p3=factor(p5);p2=(GEN)p3[1];l=lg(p2)-1;
- switch(l)
- {
- case 1: p3=gadd(gmul(gmul(p1[1],p1[2]),p1[3]),gmul(gmul(p1[4],p1[5]),p1[6]));
- p4=gsub(polx[0],p3);
- for(i=1;i<=3;i++)
- for(j=4;j<=6;j++)
- {
- p2=transroot(p1,i,j);
- p3=gadd(gmul(gmul(p2[1],p2[2]),p2[3]),gmul(gmul(p2[4],p2[5]),p2[6]));
- p4=gmul(p4,gsub(polx[0],p3));
- }
- p5=grndtoi(greal(p4),&e1);
- e1=max(e1,gexpo(gimag(p4)));
- if(e1<= 10)
- {
- p6=ggcd(p5,deriv(p5,0));
- f=(typ(p6)==10)&&(lgef(p6)>3);
- if(f) goto tchi;
- p3=factor(p5);p2=(GEN)p3[1];l=lg(p2)-1;f=carreparfait(discsr(x));
- avma=av;y=cgetg(4,17);y[3]=un;
- if(l==1)
- {
- if(f) {y[2]=un;y[1]=lstoi(360);return y;}
- else {y[2]=lneg(un);y[1]=lstoi(720);return y;}
- }
- else
- {
- if(f) {y[2]=un;y[1]=lstoi(36);return y;}
- else {y[2]=lneg(un);y[1]=lstoi(72);return y;}
- }
- }
- else prec=(prec<<2)-2;
- break;
-
- case 2: l2=lgef(p2[1])-3;if(l2>3) l2=6-l2;
- switch(l2)
- {
- case 1: f=carreparfait(discsr(x));avma=av;y=cgetg(4,17);y[3]=un;
- if(f) {y[2]=un;y[1]=lstoi(60);return y;}
- else {y[2]=lneg(un);y[1]=lstoi(120);return y;}
- case 2: f=carreparfait(discsr(x));
- if(f) {avma=av;y=cgetg(4,17);y[3]=y[2]=un;y[1]=lstoi(24);return y;}
- else
- {
- p3=(lgef(p2[1])==5)?(GEN)p2[2]:(GEN)p2[1];
- f=carreparfait(discsr(p3));avma=av;y=cgetg(4,17);y[2]=lneg(gun);
- if(f) {y[1]=lstoi(24);y[3]=deux;return y;}
- else {y[1]=lstoi(48);y[3]=un;return y;}
- }
- case 3: f=carreparfait(discsr(p2[1]))||carreparfait(discsr(p2[2]));
- avma=av;y=cgetg(4,17);y[3]=un;y[2]=lneg(gun);y[1]=lstoi(f?18:36);
- return y;
- }
- case 3: for(l2=1;l2<=3;l2++) if(lgef(p2[l2])>=6) p3=(GEN)p2[l2];
- if(lgef(p3)==6)
- {
- f=carreparfait(discsr(p3));avma=av;y=cgetg(4,17);y[2]=lneg(gun);
- y[3]=un;y[1]=f?lstoi(6):lstoi(12);return y;
- }
- else
- {
- f=carreparfait(discsr(x));avma=av;y=cgetg(4,17);y[3]=un;
- if(f) {y[2]=un;y[1]=lstoi(12);return y;}
- else {y[2]=lneg(un);y[1]=lstoi(24);return y;}
- }
- case 4: avma=av;y=cgetg(4,17);y[1]=lstoi(6);y[2]=lneg(gun);
- y[3]=deux;return y;
- default: err(galbug3);
- }
- }
- while(e1> -10);
-
- case 7:
- do
- {
- p1=roots(x,prec);p4=gun;
- for(i=1;i<=5;i++)
- for(j=i+1;j<=6;j++)
- for(k=j+1;k<=7;k++)
- p4=gmul(p4,gsub(polx[0],gadd(gadd(p1[i],p1[j]),p1[k])));
- p5=grndtoi(greal(p4),&e);e=max(e,gexpo(gimag(p4)));
- if(e> -10) prec=(prec<<2)-2;
- }
- while(e> -10);
- p6=ggcd(p5,deriv(p5,0));
- f=(typ(p6)==10)&&(lgef(p6)>3);
- if(f) goto tchi;
- p1=factor(p5);p2=(GEN)p1[1];l=lg(p2)-1;
- switch(l)
- {
- case 1: f=carreparfait(discsr(x));avma=av;y=cgetg(4,17);y[3]=un;
- if(f) {y[2]=un;y[1]=lstoi(2520);return y;}
- else {y[2]=lneg(gun);y[1]=lstoi(5040);return y;}
- case 2: f=lgef(p2[1])-3;avma=av;y=cgetg(4,17);y[3]=un;
- if((f==7)||(f==28)) {y[2]=un;y[1]=lstoi(168);return y;}
- else {y[2]=lneg(gun);y[1]=lstoi(42);return y;}
- case 3: avma=av;y=cgetg(4,17);y[3]=y[2]=un;y[1]=lstoi(21);return y;
- case 4: avma=av;y=cgetg(4,17);y[3]=un;y[2]=lneg(gun);y[1]=lstoi(14);return y;
- case 5: avma=av;y=cgetg(4,17);y[3]=y[2]=un;y[1]=lstoi(7);return y;
- default: err(galbug2);
- }
- }
- tchi:
- avma=av1;x=tchirnhausen(x1);
- }
- }
-
- GEN initalg(x,prec)
- GEN x;
- long prec;
- {
- GEN y,p1,p2,p3,p4,p5,p6,p7,fieldd,polr,ptrace,dx,adx,polmax,s,a,dxn,adxn,sn,phimax,ind;
- long tx=typ(x),n=lgef(x)-3,i,j,av=avma,av1,v,k,r1,imax,numb,flc,tetpil;
-
- if((tx!=10)||(n<=0)) err(poltyper);
- p1=content(x);p1=gcmp1(p1) ? x: gdiv(x,content(x));
- for(k=2;k<=n+2;k++) if(typ(p1[k])!=1) err(impl,"general algebraic extension");
- if(gcmp1(p2=(GEN)p1[n+2])) x=p1;
- else
- {
- x=cgetg(n+3,10);x[1]=p1[1];x[n+2]=un;x[n+1]=p1[n+1];p3=p2;
- for(k=n;k>=2;k--)
- {
- x[k]=lmulii(p3,p1[k]);
- if(k>2) p3=mulii(p2,p3);
- }
- }
- p1=factor(x);
- if(lgef(coeff(p1,1,1))!=n+3) err(redper1);
- r1=sturm(x);p4=allbase(x,0,&fieldd);
- if(r1<n)
- {
- polr=roots(x,prec);p3=cgetg(n+1,19);
- for(i=1;i<=n;i++)
- {
- p1=cgetg(n+1,18);p3[i]=(long)p1;
- for(j=1;j<=n;j++)
- p1[j]=lsubst(p4[i],varn(p4[i]),polr[j]);
- }
- p2=greal(gmul(gconj(gtrans(p3)),p3));
- }
- else
- {
- ptrace=cgetg(n+1,17);ptrace[1]=lstoi(n);
- for(k=1;k<n;k++)
- {
- p3=gmulsg(k,x[n-k+2]);
- for(i=1;i<k;i++) p3=gadd(p3,gmul(x[n-i+2],ptrace[k-i+1]));
- ptrace[k+1]=lneg(p3);
- }
- p2=cgetg(n+1,19);
- for(i=1;i<=n;i++)
- {
- p1=cgetg(n+1,18);p2[i]=(long)p1;
- for(j=1;j<i;j++) p1[j]=lcopy(coeff(p2,i,j));
- for(j=i;j<=n;j++)
- {
- p5=gres(gmul(p4[i],p4[j]),x);p6=gzero;
- for(k=0;k<=lgef(p5)-3;k++) p6=gadd(p6,gmul(p5[k+2],ptrace[k+1]));
- p1[j]=(long)p6;
- }
- }
- }
- p1=lllgram(p2,prec);v=varn(x);dx=discsr(x);adx=absi(dx);polmax=x;imax=0;
- if(r1<n) for(s=gzero,i=1;i<=n;i++) s=gadd(s,gnorm(polr[i]));
- else s=gsub(gsqr(x[n+1]),gmul2n(x[n],1));
- a=cgetg(n+1,18);for(i=1;i<=n;i++) a[i]=lmul(p4,p1[i]);
- for(numb=0,i=1;i<=n;i++)
- {
- av1=avma;p3=gmodulcp(a[i],x);p7=content(p3[2]);
- if(gcmp1(p7)) p3=caract(p3,v);
- else
- {
- p3=caract(gdiv(p3,p7),v);
- p3=gmul(gpuigs(p7,lgef(p3)-3),gsubst(p3,v,gdiv(polx[v],p7)));
- }
- p5=ggcd(deriv(p3,v),p3);
- if(lgef(p5)==3)
- {
- dxn=discsr(p3);adxn=absi(dxn);flc=gcmp(adxn,adx);numb++;
- if(flc<=0)
- {
- if(r1<n) for(sn=gzero,j=1;j<=n;j++) sn=gadd(sn,gnorm(poleval(a[i],polr[j])));
- else sn=gsub(gsqr(p3[n+1]),gmul2n(p3[n],1));
- if(flc<0) {dx=dxn;adx=adxn;s=sn;polmax=p3;imax=i;}
- else
- {
- flc=gcmp(sn,s);
- if((flc<0)||((!flc)&&(gpolcomp(p3,polmax)<0)))
- {dx=dxn;adx=adxn;s=sn;polmax=p3;imax=i;}
- }
- }
- }
- }
- if(!numb) err(polreder1);
- phimax=imax?(GEN)a[imax]:polx[v];
- j=n+1;
- while((j>=2)&&(!signe(polmax[j]))) j-=2;
- if((j>=2)&&(signe(polmax[j])>0))
- {
- for(;j>=2;j-=2) setsigne(polmax[j],-signe(polmax[j]));
- phimax=gneg(phimax);
- }
- p2=lift(gsubst(p4,v,polymodrecip(gmodulcp(phimax,x))));
- p3=cgetg(n+1,19);
- for(j=1;j<=n;j++)
- {
- p4=cgetg(n+1,18);p3[j]=(long)p4;
- for(i=1;i<=n;i++) p4[i]=(long)truecoeff(p2[j],i-1);
- }
- p4=denom(p3);
- p2=gdiv(hnf(gmul(p4,p3)),p4);p3=cgetg(n+1,17);
- for(j=1;j<=n;j++)
- {
- p1=gzero;for(i=n;i;i--) p1=gadd(coeff(p2,i,j),gmul(p1,polx[v]));
- p3[j]=(long)p1;
- }
- if(!carrecomplet(divii(dx,fieldd),&ind)) err(initalgbug1);
- tetpil=avma;y=cgetg(8,17);y[1]=lcopy(polmax);p1=cgetg(3,17);
- p1[1]=lstoi(r1);p1[2]=lstoi((n-r1)>>1);y[2]=(long)p1;
- y[3]=lcopy(fieldd);y[4]=lcopy(ind);y[5]=lcopy(s);
- if(r1<n)
- {
- p1=cgetg(n+1,18);
- for(i=1;i<=n;i++) p1[i]=(long)poleval(phimax,polr[i]);
- y[6]=(long)p1;
- }
- else y[6]=lreal(roots(polmax,prec));
- y[7]=lcopy(p3);
- return gerepile(av,tetpil,y);
- }
-
-
- GEN galoisconj2(x,prec)
- GEN x;
- long prec;
- {
- long av=avma,tetpil,i,j,n,v;
- GEN y,w,polr,p1,p2,b,di;
-
- if(typ(x)!=10) err(galer1);
- n=lgef(x)-3;if(n<=0) err(galer1);
- v=varn(x);p1=factor(x);
- if((lg(p1[1])>2)||(!gcmp1(coeff(p1,1,2)))) err(galcer1);
- polr=roots(x,prec);p1=(GEN)polr[1];b=allbase(x,0,&di);
- w=cgetg(n+1,17);w[1]=un;for(i=2;i<=n;i++) w[i]=(long)gsubst(b[i],v,p1);
- y=cgetg(n+1,17);y[1]=(long)polx[v];
- for(i=2;i<=n;i++)
- {
- p1=lindep(concat(w,polr[i]),prec);
- if(gcmp0(p1[n+1])) y[i]=zero;
- else
- {
- p2=gzero;
- for(j=1;j<=n;j++) p2=gadd(p2,gmul(p1[j],b[j]));
- p2=gdiv(p2,gneg(p1[n+1]));
- if(gcmp0(gmod(gsubst(x,v,p2),x))) y[i]=(long)p2;else y[i]=zero;
- }
- }
- tetpil=avma;return gerepile(av,tetpil,gcopy(y));
- }
-
- GEN galoisconj(x,prec)
- GEN x;
- long prec;
- {
- long av=avma,tetpil,i,j,n,v;
- GEN y,w,polr,p1,p2;
-
- if(typ(x)!=10) err(galer1);
- n=lgef(x)-3;if(n<=0) err(galer1);
- v=varn(x);p1=factor(x);
- if((lg(p1[1])>2)||(!gcmp1(coeff(p1,1,2)))) err(galcer1);
- polr=roots(x,prec);p1=(GEN)polr[1];
- w=cgetg(n+1,17);w[1]=un;for(i=2;i<=n;i++) w[i]=lmul(p1,w[i-1]);
- y=cgetg(n+1,17);y[1]=(long)polx[v];
- for(i=2;i<=n;i++)
- {
- p1=lindep(concat(w,polr[i]),prec);
- if(gcmp0(p1[n+1])) y[i]=zero;
- else
- {
- p2=gzero;
- for(j=n;j;j--) p2=gadd(p1[j],gmul(p2,polx[v]));
- p2=gdiv(p2,gneg(p1[n+1]));
- if(gcmp0(gmod(gsubst(x,v,p2),x))) y[i]=(long)p2;else y[i]=zero;
- }
- }
- tetpil=avma;return gerepile(av,tetpil,gcopy(y));
- }
-
-
-
-
-
-
-
-
-
-
-