home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Professional / OS2PRO194.ISO / os2 / progs / pari / pari_137 / src / base.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-05-20  |  35.3 KB  |  1,433 lines
  1.  /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  2. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  3. /*                                                                 */
  4. /*                       BASE D'ENTIERS                            */
  5. /*                                                                 */
  6. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  7. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  8.  
  9. # include "genpari.h"
  10. GEN rquot(),ordmax(),rtran(),mtran(),matinv();
  11. void rowred();
  12.  
  13. GEN allbase(x,code,y)
  14.  
  15.      GEN x,*y;
  16.      long code;
  17.  
  18. /*******************************************************************
  19.   Entree:     x polynome unitaire a coefficients dans Z de deg n
  20.         definissant un corps de nombres K=Q(theta);
  21.               code 0, 1 ou (long)p selon que l'on veut base, smallbase
  22.             ou factoredbase.
  23.           y pointeur sur un GEN destine a recevoir
  24.         le discriminant du corps K.
  25.   Sortie:    retourne 1) un vecteur (horizontal) a n composantes, 
  26.             de polynomes a coeff dans Q (de deg 0,1...n-1)
  27.         constituant une base de l'anneau des entiers de K.
  28.             2) le discriminant de K (dans *y).
  29.         Rem: le denominateur commun des coef. est dans da.
  30. *******************************************************************/
  31.  
  32. {
  33.   GEN p,a,at,bt,b,da,db,q;
  34.   long av=avma,tetpil,n,h,j,je,k,r,s,t,pro,v;
  35.  
  36.   if(typ(x)!=10) err(allbaser1);
  37.   n=lgef(x)-3;if(n<=0) err(allbaser1);
  38.   v=varn(x);*y=discsr(x);
  39.   switch(code)
  40.     {
  41.     case 0: p=auxdecomp(absi(*y),1);h=lg(p[1])-1;break; /* base */
  42.     case 1: p=auxdecomp(absi(*y),0);h=lg(p[1])-1;break; /* smallbase */
  43.     default: p=(GEN)code;
  44.       if((typ(p)!=19)||(lg(p)!=3)) err(factoreder1); /* factoredbase */
  45.       h=lg(p[1])-1;
  46.       q=gun;for(je=1;je<=h;je++) q=gmul(q,gpui(coeff(p,je,1),coeff(p,je,2)));
  47.       if(gcmp(absi(q),absi(*y))) err(factoreder2);
  48.     }
  49.   a=idmat(n);da=gun;
  50.   for(je=1;je<=h;je++)
  51.     {
  52.       if(gcmpgs(coeff(p,je,2),1)>0)
  53.     {
  54.       b=ordmax(x,coeff(p,je,1),coeff(p,je,2),&db);
  55.       a=gmul(db,a);b=gmul(da,b);
  56.       da=mulii(db,da);db=da;
  57.       at=gtrans(a);bt=gtrans(b);
  58.       for(r=n-1;r>=0;r--)
  59.         for(s=r;s>=0;s--)
  60.           while(signe(coef1(bt,s,r)))
  61.         {
  62.           q=rquot(coef1(at,s,s),coef1(bt,s,r));
  63.           at[s+1]=(long)rtran(at[s+1],bt[r+1],q);
  64.           for(t=s-1;t>=0;t--)
  65.             {
  66.               q=rquot(coef1(at,t,s),coef1(at,t,t));
  67.               at[s+1]=(long)rtran(at[s+1],at[t+1],q);
  68.             }
  69.           pro=at[s+1];at[s+1]=bt[r+1];bt[r+1]=pro;
  70.         }
  71.       for(j=n-1;j>=0;j--)
  72.         {
  73.           for(k=0;k<j;k++)
  74.         {
  75.           while(signe(coef1(at,j,k)))
  76.             {
  77.               q=rquot(coef1(at,j,j),coef1(at,j,k));
  78.               at[j+1]=(long)rtran(at[j+1],at[k+1],q);
  79.               pro=at[j+1];at[j+1]=at[k+1];at[k+1]=pro;
  80.             }
  81.         }
  82.           if(signe(coef1(at,j,j))<0)
  83.         for(k=0;k<=j;k++) coef1(at,k,j)=lnegi(coef1(at,k,j));
  84.           for(k=j+1;k<n;k++)
  85.         {
  86.           q=rquot(coef1(at,j,k),coef1(at,j,j));
  87.           at[k+1]=(long)rtran(at[k+1],at[j+1],q);
  88.         }
  89.         }
  90.       for(j=1;j<n;j++)
  91.         if(!cmpii(coef1(at,j,j),coef1(at,j-1,j-1)))
  92.           {
  93.         coef1(at,0,j)=zero;
  94.         for(k=1;k<=j;k++)
  95.           coef1(at,k,j)=coef1(at,k-1,j-1);
  96.           }
  97.       a=gtrans(at);
  98.     }
  99.     }
  100.   for(j=1;j<=n;j++)
  101.     {
  102.       *y=divii(mulii(coeff(a,j,j),*y),da);
  103.       *y=divii(mulii(coeff(a,j,j),*y),da);
  104.     }
  105.   tetpil=avma;*y=gcopy(*y);at=cgetg(n+1,17);
  106.   for(j=1;j<=n;j++)
  107.     {
  108.       q=cgetg(j+2,10);q[1]=0x1000002+j+(v<<16);at[j]=(long)q;
  109.       for(k=2;k<=j+1;k++) q[k]=ldiv(coeff(a,j,k-1),da);
  110.     }
  111.   pro=lpile(av,tetpil,0)>>2;at+=pro;(*y)+=pro;
  112.   return at;
  113. }
  114.  
  115. GEN base(x,y)
  116.      GEN x,*y;
  117. {
  118.   return allbase(x,0,y);
  119. }
  120.  
  121. GEN smallbase(x,y)
  122.      GEN x,*y;
  123. {
  124.   return allbase(x,1,y);
  125. }
  126.  
  127. GEN factoredbase(x,p,y)
  128.      GEN x,p,*y;
  129. {
  130.   return allbase(x,(long)p,y);
  131. }
  132.  
  133. GEN discf(x)
  134.      GEN x;
  135. {
  136.   GEN y;
  137.   long av,tetpil;
  138.  
  139.   av=avma;allbase(x,0,&y);tetpil=avma;
  140.   return gerepile(av,tetpil,gcopy(y));
  141. }
  142.  
  143. GEN smalldiscf(x)
  144.      GEN x;
  145. {
  146.   GEN y;
  147.   long av,tetpil;
  148.  
  149.   av=avma;allbase(x,1,&y);tetpil=avma;
  150.   return gerepile(av,tetpil,gcopy(y));
  151. }
  152.  
  153. GEN factoreddiscf(x,p)
  154.      GEN x,p;
  155. {
  156.   GEN y;
  157.   long av,tetpil;
  158.  
  159.   av=avma;allbase(x,(long)p,&y);tetpil=avma;
  160.   return gerepile(av,tetpil,gcopy(y));
  161. }
  162.  
  163. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  164. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  165. /*                                                                 */
  166. /*   Quotient et Reste normalises   ( -1/2 < r = x-q*y <= 1/2 )    */
  167. /*                                                                 */
  168. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  169. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  170.  
  171. GEN rquot(x,y)
  172.  
  173.      GEN x,y;
  174.  
  175. {
  176.   GEN u,v,w,p;
  177.   long av,av1;
  178.  
  179.   av=avma;
  180.   u=absi(y);v=shifti(x,1);w=shifti(y,1);
  181.   if ( cmpii(u,v)>0) p=subii(v,u);
  182.   else p=addsi(-1,addii(u,v));
  183.   av1=avma;
  184.   return(gerepile(av,av1,divii(p,w)));
  185. }
  186.  
  187. GEN rrmdr(x,y)
  188.  
  189.      GEN x,y;
  190.  
  191. {
  192.   GEN p;
  193.   long av,av1;
  194.  
  195.   av=avma;
  196.   p=mulii(rquot(x,y),y);
  197.   av1=avma;
  198.   return(gerepile(av,av1,subii(x,p)));
  199. }
  200.  
  201. GEN rinv(x,y)
  202.  
  203.      GEN x,y;
  204.  
  205. {
  206.   GEN a,c,q,r,t;
  207.   long av,av1;
  208.  
  209.   av=avma;
  210.   a=gun;c=gzero;
  211.   while( signe(y))
  212.     {
  213.       q=rquot(x,y);
  214.       r=subii(a,mulii(q,c));a=c;c=r;
  215.       t=subii(x,mulii(q,y));x=y;y=t;
  216.     }
  217.   av1=avma;
  218.   if (signe(x)<0) a=negi(a);
  219.   if (signe(c)) { av1=avma; a=rrmdr(a,c); }
  220.   return( gerepile(av,av1,a));
  221. }
  222.  
  223. GEN rgcd(x,y)
  224.  
  225.      GEN x,y;
  226.  
  227. {
  228.   GEN z;
  229.   long av,av1;
  230.  
  231.   av=avma;
  232.   while(signe(y))
  233.     {
  234.       z=rrmdr(x,y);x=y;y=z;
  235.     }
  236.   av1=avma;
  237.   return(gerepile(av,av1,absi(x)));
  238. }
  239.     
  240.  
  241. GEN rlcm(x,y)
  242.  
  243.      GEN x,y;
  244.  
  245. {
  246.   GEN d,z;
  247.   long av,av1;
  248.  
  249.   av=avma;
  250.   z=mulii(x,y);d=rgcd(x,y);
  251.   av1=avma;
  252.   return(gerepile(av,av1,divii(z,d)));
  253. }
  254.  
  255. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  256. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  257. /*                                                                 */
  258. /*           Matrice compagnon du polynome unitaire x              */
  259. /*                                                                 */
  260. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  261. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  262.  
  263. GEN companion(x)
  264.  
  265.      GEN     x;
  266.  
  267. {
  268.   long    i,j,l;
  269.   GEN     y;
  270.   
  271.   l=lgef(x)-2;y=cgetg(l,19);
  272.   for(i=1;i<l;i++) y[i]=lgetg(l,18);
  273.   for(i=0;i<l-2;i++)
  274.     for(j=0;j<l-1;j++) coef1(y,i,j)=((i+1)==j) ? un : zero;
  275.   for(j=0;j<l-1;j++) coef1(y,l-2,j)=lneg(x[j+2]);
  276.   return(y);
  277. }
  278.  
  279.  
  280.  
  281. GEN ordmax(f,p,e,ptdelta)
  282.  
  283. GEN f,p,e;
  284. GEN *ptdelta;
  285.  
  286. {
  287.   GEN m,hh,pp,dd,ppdd,index,q,r,s,b,c,t,jp,v,delta;
  288.   GEN cf[20],w[20],a;
  289.   long h,i,j,k,sp,epsilon,n=lgef(f)-3,av=avma,tetpil,dec;
  290.  
  291.   a=cgetg(n*n+1,19);
  292.   for(j=1;j<=n*n;j++)
  293.   {
  294.     a[j]=lgetg(n+1,18);
  295.     for(k=1;k<=n;k++) coeff(a,k,j)=zero;
  296.   }
  297.   v=cgetg(n+1,18);
  298.   cf[0]=idmat(n);
  299.   cf[1]=companion(f);
  300.   for(j=2;j<n;j++) cf[j]=gmul(cf[1],cf[j-1]);
  301.   delta=gun; epsilon=itos(e);
  302.   m=idmat(n);
  303.  
  304.   do
  305.     {
  306.       pp=mulii(p,p);
  307.       dd=mulii(delta,delta);
  308.       ppdd=mulii(dd,pp);
  309.       b=matinv(m,delta);
  310.       for(i=0;i<n;i++)
  311.     {
  312.       t=gscalsmat(0,n); /* t <--- matrice nulle d'ordre n */
  313.       for(h=0;h<n;h++)
  314.         for(j=0;j<n;j++)
  315.           for(k=0;k<n;k++)
  316.         coef1(t,j,k)=(long)rrmdr(addii(coef1(t,j,k),mulii(coef1(m,i,h),coef1(cf[h],j,k))),ppdd);
  317.       c=gmul(t,b);
  318.       w[i]=gmul(m,c);
  319.       for(j=0;j<n;j++)
  320.         for(k=0;k<n;k++)
  321.           coef1(w[i],j,k)=(long)rrmdr(divii(coef1(w[i],j,k),dd),pp);
  322.     }
  323.       if(cmpis(p,n)>0)
  324.     {
  325.       for(i=0;i<n;i++)
  326.         for(j=0;j<n;j++)
  327.           {
  328.         coeff(t,i+1,j+1)=zero;
  329.         for(k=0;k<n;k++)
  330.           for(h=0;h<n;h++)
  331.             {
  332.               r=modii(coef1(w[i],k,h),p);
  333.               s=modii(coef1(w[j],h,k),p);
  334.               coef1(t,i,j)=lmodii(addii(coef1(t,i,j),mulii(r,s)),p);
  335.             }
  336.           }
  337.     }
  338.       else
  339.     {
  340.       for(j=0;j<n;j++)
  341.         {
  342.           for(i=0;i<n;i++)
  343.         coef1(b,i,j)=(i==0)? un:zero;
  344. /* ici la boucle en k calcule la puissance p mod p de w[j] */
  345.           sp=itos(p);
  346.           for(k=0;k<sp;k++)
  347.         {
  348.           for(i=0;i<n;i++)
  349.             {
  350.               v[i+1]=zero;
  351.               for(h=0;h<n;h++)
  352.             v[i+1]=lmodii(addii(v[i+1],mulii(coef1(b,h,j),coef1(w[j],h,i))),p);
  353.             }
  354.           for(i=0;i<n;i++) coef1(b,i,j)=v[i+1];
  355.         }
  356.         }
  357.       q=p;t=b;
  358.       while(cmpis(q,n)<0)
  359.         {
  360.           q=mulii(q,p);
  361.           t=gmul(b,t);
  362.         }
  363.     }
  364.       for(i=0;i<n;i++)
  365.     for(j=0;j<n;j++)
  366.       {
  367.         coef1(a,j,i)=(i==j)? (long)p:zero;
  368.         coef1(a,j,n+i)=lmodii(coef1(t,i,j),p);
  369.       }
  370.       rowred(a,2*n-1,pp);
  371.       for(i=0;i<n;i++)
  372.     for(j=0;j<n;j++)
  373.       coef1(b,j,i)=coef1(a,j,i);
  374.       jp=matinv(b,p);
  375.       for(k=0;k<n;k++)
  376.     {
  377.       t=gmul(jp,w[k]);
  378.       t=gmul(t,b);
  379.       for(i=0;i<n;i++)
  380.         for(j=0;j<n;j++)
  381.           coef1(t,i,j)=ldivii(coef1(t,i,j),p);
  382.       h=0;
  383.       for(i=0;i<n;i++)
  384.         for(j=0;j<n;j++)
  385.           {
  386.         coef1(a,k,h)=coef1(t,i,j);
  387.         h++;
  388.           }
  389.     }
  390.       rowred(a,n*n-1,pp);
  391.       index=gun;
  392.       for(i=0;i<n;i++)
  393.     index=mulii(index,coef1(a,i,i));
  394.       if (cmpsi(1,index))
  395.     {
  396.       delta=mulii(index,delta);
  397.       for(i=0;i<n;i++)
  398.         for(j=0;j<n;j++)
  399.           coef1(c,i,j)=coef1(a,i,j);
  400.       b=matinv(c,index);
  401.       m=gmul(b,m);
  402.       hh=delta;
  403.       for(i=0;i<n;i++)
  404.         for(j=0;j<n;j++)
  405.           hh=rgcd(coef1(m,i,j),hh);
  406.       if(cmpis(hh,1)>1)
  407.         {
  408.           delta=divii(delta,hh);
  409.           for(i=0;i<n;i++)
  410.         for(j=0;j<n;j++)
  411.           coef1(m,i,j)=ldivii(coef1(m,i,j),hh);
  412.         }
  413.       q=index;
  414.       while(!signe(modii(q,p)))
  415.         {
  416.           q=divii(q,p);
  417.           epsilon=epsilon-2;
  418.         }
  419.     }
  420.     }
  421.   while(!gcmp1(index) && (epsilon>=2));
  422.   tetpil=avma;delta=gcopy(delta);m=gcopy(m);
  423.   dec=lpile(av,tetpil,0)>>2;
  424.   *ptdelta=delta+dec;
  425.   return m+dec;
  426. }
  427.  
  428. void rowred(a,rlim,rmod)
  429.      GEN a,rmod;
  430.      long rlim;
  431.  
  432. {
  433.   long j,k,n,pro;
  434.   GEN q;
  435.  
  436.   n=lg(a[1])-1;
  437.   for(j=0;j<n;j++)
  438.     {
  439.       for(k=j+1;k<=rlim;k++)
  440.     while (signe(coef1(a,j,k)))
  441.       {
  442.         q=rquot(coef1(a,j,j),coef1(a,j,k));
  443.         a[j+1]=(long)mtran(a[j+1],a[k+1],q,rmod);
  444.         pro=a[j+1];a[j+1]=a[k+1];a[k+1]=pro;
  445.       }
  446.       if (signe(coef1(a,j,j))<0)
  447.     for(k=j;k<n;k++) coef1(a,k,j)=lnegi(coef1(a,k,j));
  448.       for(k=0;k<j;k++)
  449.     {
  450.       q=rquot(coef1(a,j,k),coef1(a,j,j));
  451.       a[k+1]=(long)mtran(a[k+1],a[j+1],q,rmod);
  452.     }
  453.     }
  454. }
  455.  
  456. GEN rtran(v,w,q)
  457.      GEN v,w,q ;
  458.  
  459. {
  460.   long av,tetpil;
  461.   GEN p1;
  462.  
  463.   if (signe(q))
  464.     {
  465.       av=avma;p1=gmul(q,w);tetpil=avma;
  466.       return gerepile(av,tetpil,gsub(v,p1));
  467.     }
  468.   else return v;
  469. }
  470.  
  471. GEN mtran(v,w,q,m)
  472.      GEN v,w,q,m;
  473.  
  474. {
  475.   long k;
  476.   
  477.   if (signe(q))
  478.     {
  479.       for(k=0;k<lg(v)-1;k++)
  480.     {
  481.       v[k+1]=(long)rrmdr(subii(v[k+1],modii(mulii(q,w[k+1]),m)),m);
  482.     }
  483.     }
  484.   return v;
  485. }
  486.  
  487.  
  488. GEN matinv(x,d)
  489.      GEN x,d;
  490. /*=======================================================================
  491.     Calcule d/x  ou  d est entier et x matrice triangulaire inferieure
  492.   entiere dont les coeff diagonaux divisent d ( resultat entier).
  493. ========================================================================*/
  494. {
  495.   long n,h,i,j,k,av,av1;
  496.   GEN y;
  497.  
  498.   av=avma;
  499.   y=idmat(n=lg(x)-1);
  500.   for(i=1;i<=n;i++)
  501.     coeff(y,i,i)=ldivii(d,coeff(x,i,i));
  502.   for(i=2;i<=n;i++)
  503.     for(j=i-1;j;j--)
  504.       {
  505.     for(h=zero,k=j+1;k<=i;k++)
  506.       h=ladd(h,mulii(coeff(y,i,k),coeff(x,k,j)));
  507.     coeff(y,i,j)=ldivii(negi(h),coeff(x,j,j));
  508.       }
  509.   av1=avma;
  510.   return gerepile(av,av1,gcopy(y));
  511. }
  512.  
  513. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  514. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  515. /*~                                    ~*/
  516. /*~            HERMITE NORMAL FORM REDUCTION            ~*/
  517. /*~                                    ~*/
  518. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  519. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  520.  
  521. GEN hnfold(x)
  522.      GEN x;
  523. {
  524.   long li,co,av,tetpil,i,j,j1,def,ind,lim;
  525.   GEN p1,p2,y,z,m1;
  526.  
  527.   if(typ(x)!=19) err(hnfer1);
  528.   lim=(avma+bot)>>1;
  529.   av=avma;co=lg(x);li=lg(x[1]);y=gcopy(x);
  530.   def=co;
  531.   if(li>co) err(hnfer2);
  532.   for(i=li-1;i>=1;i--)
  533.     {
  534.       def--;j=1;while((j<def)&&(!signe(coeff(y,i,j)))) j++;
  535.       while(j<def)
  536.     {
  537.       m1=absi(coeff(y,i,j));ind=j;
  538.       for(j1=j+1;j1<def;j1++)
  539.         {
  540.           p1=(GEN)coeff(y,i,j1);
  541.           if(signe(p1)&&(cmpii(absi(p1),m1)<0))
  542.         {
  543.           m1=p1;ind=j1;
  544.         }
  545.         }
  546.       p1=(GEN)coeff(y,i,def);
  547.       if(!(signe(p1)&&(cmpii(absi(p1),m1)<0)))
  548.         {
  549.           p1=(GEN)y[def];y[def]=y[ind];y[ind]=(long)p1;
  550.         }
  551.       p1=(GEN)coeff(y,i,def);if(signe(p1)<0) y[def]=lneg(y[def]);
  552.       p1=(GEN)coeff(y,i,def);
  553.       for(j=1;j<def;j++)
  554.         {
  555.           p2=gdivent(coeff(y,i,j),p1);
  556.           y[j]=lsub(y[j],gmul(p2,y[def]));
  557.         }
  558.       j=1;while((j<def)&&(!signe(coeff(y,i,j)))) j++;
  559.     }
  560.       p1=(GEN)coeff(y,i,def);
  561.       if(signe(p1)<0) {y[def]=lneg(y[def]);p1=(GEN)coeff(y,i,def);}
  562.       if(signe(p1))
  563.     {
  564.       for(j=def+1;j<co;j++)
  565.         {
  566.           p2=gdivent(coeff(y,i,j),p1);
  567.           y[j]=lsub(y[j],gmul(p2,y[def]));
  568.         }
  569.     }
  570.       if(avma<lim) {tetpil=avma;y=gerepile(av,tetpil,gcopy(y));}
  571.     }
  572.   tetpil=avma;z=cgetg(li,19);
  573.   for(j=1;j<li;j++)
  574.     {
  575.       z[j]=lcopy(y[j+co-li]);
  576.     }
  577.   for(j=1;j<=co-li;j++) if(!gcmp0(y[j])) err(hnfer3);
  578.   return gerepile(av,tetpil,z);
  579. }
  580.  
  581. GEN hnfold2(x)
  582.      GEN x;
  583. {
  584.   long li,co,av,tetpil,i,j,def,lim;
  585.   GEN p1,p2,y,z,u,v,d;
  586.  
  587.   if(typ(x)!=19) err(hnfer1);
  588.   lim=(avma+bot)>>1;
  589.   av=avma;co=lg(x);li=lg(x[1]);y=gcopy(x);
  590.   def=co;
  591.   if(li>co) err(hnfer2);
  592.   for(i=li-1;i>=1;i--)
  593.     {
  594.       def--;j=def;while(j&&(!signe(coeff(y,i,j)))) j--;
  595.       while(j>1)
  596.     {
  597.       d=bezout(coeff(y,i,j),coeff(y,i,j-1),&u,&v);
  598.       p1=gadd(gmul(u,y[j]),gmul(v,y[j-1]));
  599.       y[j]=lsub(gmul(divii(coeff(y,i,j),d),y[j-1]),gmul(divii(coeff(y,i,j-1),d),y[j]));
  600.       y[j-1]=(long)p1;
  601.       j--;while(j&&(!signe(coeff(y,i,j)))) j--;
  602.     }
  603.       if(j==1)
  604.     {
  605.       d=bezout(coeff(y,i,1),coeff(y,i,def),&u,&v);
  606.       p1=gadd(gmul(u,y[1]),gmul(v,y[def]));
  607.       y[1]=lsub(gmul(divii(coeff(y,i,1),d),y[def]),gmul(divii(coeff(y,i,def),d),y[1]));
  608.       y[def]=(long)p1;
  609.     }
  610.       p1=(GEN)coeff(y,i,def);
  611.       if(signe(p1)<0) {y[def]=lneg(y[def]);p1=(GEN)coeff(y,i,def);}
  612.       if(signe(p1))
  613.     {
  614.       for(j=def+1;j<co;j++)
  615.         {
  616.           p2=gdivent(coeff(y,i,j),p1);
  617.           y[j]=lsub(y[j],gmul(p2,y[def]));
  618.         }
  619.     }
  620.       if(avma<lim) {tetpil=avma;y=gerepile(av,tetpil,gcopy(y));}
  621.     }
  622.   tetpil=avma;z=cgetg(li,19);
  623.   for(j=1;j<li;j++)
  624.     {
  625.       z[j]=lcopy(y[j+co-li]);
  626.     }
  627.   for(j=1;j<=co-li;j++) if(!gcmp0(y[j])) err(hnfer3);
  628.   return gerepile(av,tetpil,z);
  629. }
  630.  
  631. GEN hnf(x)
  632.      GEN x;
  633. {
  634.   long li,co,av,tetpil,i,j,k,def,ldef,lim;
  635.   GEN p1,p2,y,z,u,v,d;
  636.  
  637.   if(typ(x)!=19) err(hnfer1);
  638.   lim=(avma+bot)>>1;
  639.   av=avma;co=lg(x);li=lg(x[1]);y=gcopy(x);
  640.   def=co;ldef=(li>co)?li-co+1:1;
  641.   for(i=li-1;i>=ldef;i--)
  642.     {
  643.       def--;j=def-1;while(j&&(!signe(coeff(y,i,j)))) j--;
  644.       while(j>1)
  645.     {
  646.       d=bezout(coeff(y,i,j),coeff(y,i,j-1),&u,&v);
  647.       p1=gadd(gmul(u,y[j]),gmul(v,y[j-1]));
  648.       y[j]=lsub(gmul(divii(coeff(y,i,j),d),y[j-1]),gmul(divii(coeff(y,i,j-1),d),y[j]));
  649.       y[j-1]=(long)p1;
  650.       j--;while(j&&(!signe(coeff(y,i,j)))) j--;
  651.     }
  652.       if(j==1)
  653.     {
  654.       d=bezout(coeff(y,i,1),coeff(y,i,def),&u,&v);
  655.       p1=gadd(gmul(u,y[1]),gmul(v,y[def]));
  656.       y[1]=lsub(gmul(divii(coeff(y,i,1),d),y[def]),gmul(divii(coeff(y,i,def),d),y[1]));
  657.       y[def]=(long)p1;
  658.     }
  659.       p1=(GEN)coeff(y,i,def);
  660.       if(signe(p1)<0) {y[def]=lneg(y[def]);p1=(GEN)coeff(y,i,def);}
  661.       if(signe(p1))
  662.     {
  663.       for(j=def+1;j<co;j++)
  664.         {
  665.           p2=gdivent(coeff(y,i,j),p1);
  666.           y[j]=lsub(y[j],gmul(p2,y[def]));
  667.         }
  668.     }
  669.       else def++;
  670.       if(avma<lim) {tetpil=avma;y=gerepile(av,tetpil,gcopy(y));}
  671.     }
  672.   for(i=0,j=1;j<co;j++) if(!gcmp0(y[j])) i++;
  673.   tetpil=avma;z=cgetg(i+1,19);
  674.   for(k=0,j=1;j<co;j++) if(!gcmp0(y[j])) z[++k]=lcopy(y[j]);
  675.   return gerepile(av,tetpil,z);
  676. }
  677.  
  678.  
  679. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  680. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  681. /*~                                    ~*/
  682. /*~            SMITH NORMAL FORM REDUCTION            ~*/
  683. /*~                                    ~*/
  684. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  685. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  686.  
  687. GEN smith2(x)
  688.      GEN x;
  689. {
  690.   long li,av,tetpil,i,j,k,l,lim,c,fl,n;
  691.   GEN p1,p2,p3,p4,y,z,b,u,v,d;
  692.  
  693.   if(typ(x)!=19) err(hnfer1);
  694.   lim=(avma+bot)>>1;
  695.   av=avma;n=lg(x)-1;li=lg(x[1])-1;y=gcopy(x);
  696.   if(li!=n) err(hnfer2);
  697.   for(i=n;i>=2;i--)
  698.     {
  699.       do
  700.     {
  701.       for(j=i-1;j>=1;j--)
  702.         {
  703.           p1=(GEN)coeff(y,i,j);
  704.           if(signe(p1))
  705.         {
  706.           p2=(GEN)coeff(y,i,i);
  707.           if(gegal(p1,p2)) {d=p1;u=gun;v=gzero;p3=gun;p4=gun;}
  708.           else if(!signe(addii(p1,p2))) 
  709.             {
  710.               d=absi(p1);u=(signe(p2)>0)?gun:gneg(un);
  711.               v=gzero;p3=u;p4=gneg(u);
  712.             }
  713.           else {d=bezout(p2,p1,&u,&v);p3=divii(p2,d);p4=divii(p1,d);}
  714.           for(k=1;k<=i;k++)
  715.             {
  716.               b=addii(mulii(u,coeff(y,k,i)),mulii(v,coeff(y,k,j)));
  717.               coeff(y,k,j)=lsubii(mulii(p3,coeff(y,k,j)),mulii(p4,coeff(y,k,i)));
  718.               coeff(y,k,i)=(long)b;
  719.             }
  720.         }
  721.         }
  722.       c=0;
  723.       for(j=i-1;j>=1;j--)
  724.         {
  725.           p1=(GEN)coeff(y,j,i);
  726.           if(signe(p1))
  727.         {
  728.           p2=(GEN)coeff(y,i,i);
  729.           if(gegal(p1,p2)) {d=p1;u=gun;v=gzero;p3=gun;p4=gun;}
  730.           else if(!signe(addii(p1,p2))) 
  731.             {
  732.               d=absi(p1);u=(signe(p2)>0)?gun:gneg(un);
  733.               v=gzero;p3=u;p4=gneg(u);
  734.             }
  735.           else {d=bezout(p2,p1,&u,&v);p3=divii(p2,d);p4=divii(p1,d);}
  736.           for(k=1;k<=i;k++)
  737.             {
  738.               b=addii(mulii(u,coeff(y,i,k)),mulii(v,coeff(y,j,k)));
  739.               coeff(y,j,k)=lsubii(mulii(p3,coeff(y,j,k)),mulii(p4,coeff(y,i,k)));
  740.               coeff(y,i,k)=(long)b;
  741.             }
  742.           c++;
  743.         }
  744.         }
  745.       if(!c)
  746.         {
  747.           b=(GEN)coeff(y,i,i);fl=1;
  748.           if(signe(b))
  749.         {
  750.           for(k=1;(k<i)&&fl;k++)
  751.             for(l=1;(l<i)&&fl;l++)
  752.               fl= !signe(modii(coeff(y,k,l),b));
  753.         }
  754.           if(!fl) {l--;y[i]=ladd(y[i],y[l]);}
  755.         }
  756.       if(avma<lim) {tetpil=avma;y=gerepile(av,tetpil,gcopy(y));}
  757.     }
  758.       while(c||(!fl));
  759.     }
  760.   tetpil=avma;z=cgetg(n+1,17);
  761.   for(k=1;k<=n;k++) z[k]=labs(coeff(y,k,k));
  762.   return gerepile(av,tetpil,z);
  763. }
  764.  
  765. GEN smith(x)
  766.      GEN x;
  767. {
  768.   long li,av,tetpil,i,j,k,l,lim,c,fl,n;
  769.   GEN p1,p2,p3,p4,y,z,b,u,v,d;
  770.  
  771.   if(typ(x)!=19) err(hnfer1);
  772.   lim=(avma+bot)>>1;
  773.   av=avma;n=lg(x)-1;if(!n) return cgetg(1,17);
  774.   li=lg(x[1])-1;y=gcopy(x);
  775.   if(li!=n) err(hnfer2);
  776.   for(i=n;i>=2;i--)
  777.     {
  778.       do
  779.     {
  780.       c=0;
  781.       for(j=i-1;j>=1;j--)
  782.         {
  783.           p1=(GEN)coeff(y,i,j);
  784.           if(signe(p1))
  785.         {
  786.           p2=(GEN)coeff(y,i,i);
  787.           if(gegal(p1,p2)) {d=p1;u=gun;v=gzero;p3=gun;p4=gun;}
  788.           else if(!signe(addii(p1,p2))) 
  789.             {
  790.               d=absi(p1);u=(signe(p2)>0)?gun:gneg(un);
  791.               v=gzero;p3=u;p4=gneg(u);
  792.             }
  793.           else {d=bezout(p2,p1,&u,&v);p3=divii(p2,d);p4=divii(p1,d);}
  794.           for(k=1;k<=i;k++)
  795.             {
  796.               b=addii(mulii(u,coeff(y,k,i)),mulii(v,coeff(y,k,j)));
  797.               coeff(y,k,j)=lsubii(mulii(p3,coeff(y,k,j)),mulii(p4,coeff(y,k,i)));
  798.               coeff(y,k,i)=(long)b;
  799.             }
  800.           c++;
  801.         }
  802.         }
  803.       for(j=i-1;j>=1;j--)
  804.         {
  805.           p1=(GEN)coeff(y,j,i);
  806.           if(signe(p1))
  807.         {
  808.           p2=(GEN)coeff(y,i,i);
  809.           if(gegal(p1,p2)) {d=p1;u=gun;v=gzero;p3=gun;p4=gun;}
  810.           else if(!signe(addii(p1,p2))) 
  811.             {
  812.               d=absi(p1);u=(signe(p2)>0)?gun:gneg(un);
  813.               v=gzero;p3=u;p4=gneg(u);
  814.             }
  815.           else {d=bezout(p2,p1,&u,&v);p3=divii(p2,d);p4=divii(p1,d);}
  816.           for(k=1;k<=i;k++)
  817.             {
  818.               b=addii(mulii(u,coeff(y,i,k)),mulii(v,coeff(y,j,k)));
  819.               coeff(y,j,k)=lsubii(mulii(p3,coeff(y,j,k)),mulii(p4,coeff(y,i,k)));
  820.               coeff(y,i,k)=(long)b;
  821.             }
  822.         }
  823.         }
  824.       if(!c)
  825.         {
  826.           b=(GEN)coeff(y,i,i);fl=1;
  827.           if(signe(b))
  828.         {
  829.           for(k=1;(k<i)&&fl;k++)
  830.             for(l=1;(l<i)&&fl;l++)
  831.               fl= !signe(modii(coeff(y,k,l),b));
  832.         }
  833.           if(!fl)
  834.         {
  835.           k--;
  836.           for(l=1;l<=i;l++)
  837.             coeff(y,i,l)=laddii(coeff(y,i,l),coeff(y,k,l));
  838.         }
  839.         }
  840.       if(avma<lim) {tetpil=avma;y=gerepile(av,tetpil,gcopy(y));}
  841.     }
  842.       while(c||(!fl));
  843.     }
  844.   tetpil=avma;z=cgetg(n+1,17);
  845.   for(j=0,k=1;k<=n;k++) if(!signe(coeff(y,k,k))) z[++j]=zero;
  846.   for(k=1;k<=n;k++) if(signe(coeff(y,k,k))) z[++j]=labs(coeff(y,k,k));
  847.   return gerepile(av,tetpil,z);
  848. }
  849.  
  850.  
  851. GEN oldsmith(x)
  852.      GEN x;
  853. {
  854.   long li,av,tetpil,i,j,k,l,lim,c,fl,n;
  855.   GEN p1,p2,p3,p4,y,z,b,u,v,d;
  856.  
  857.   if(typ(x)!=19) err(hnfer1);
  858.   lim=(avma+bot)>>1;
  859.   av=avma;n=lg(x)-1;li=lg(x[1])-1;y=gcopy(x);
  860.   if(li!=n) err(hnfer2);
  861.   for(i=n;i>=2;i--)
  862.     {
  863.       do
  864.     {
  865.       c=0;
  866.       for(j=i-1;j>=1;j--)
  867.         {
  868.           p1=(GEN)coeff(y,i,j);
  869.           if(signe(p1))
  870.         {
  871.           p2=(GEN)coeff(y,i,i);
  872.           if(gegal(p1,p2)) {d=p1;u=gun;v=gzero;p3=gun;p4=gun;}
  873.           else {d=bezout(p2,p1,&u,&v);p3=divii(p2,d);p4=divii(p1,d);}
  874.           for(k=1;k<=i;k++)
  875.             {
  876.               b=addii(mulii(u,coeff(y,k,i)),mulii(v,coeff(y,k,j)));
  877.               coeff(y,k,j)=lsubii(mulii(p3,coeff(y,k,j)),mulii(p4,coeff(y,k,i)));
  878.               coeff(y,k,i)=(long)b;
  879.             }
  880.           c++;
  881.         }
  882.         }
  883.       for(j=i-1;j>=1;j--)
  884.         {
  885.           p1=(GEN)coeff(y,j,i);
  886.           if(signe(p1))
  887.         {
  888.           p2=(GEN)coeff(y,i,i);
  889.           if(gegal(p1,p2)) {d=p1;u=gun;v=gzero;p3=gun;p4=gun;}
  890.           else {d=bezout(p2,p1,&u,&v);p3=divii(p2,d);p4=divii(p1,d);}
  891.           for(k=1;k<=i;k++)
  892.             {
  893.               b=addii(mulii(u,coeff(y,i,k)),mulii(v,coeff(y,j,k)));
  894.               coeff(y,j,k)=lsubii(mulii(p3,coeff(y,j,k)),mulii(p4,coeff(y,i,k)));
  895.               coeff(y,i,k)=(long)b;
  896.             }
  897.         }
  898.         }
  899.       if(!c)
  900.         {
  901.           b=(GEN)coeff(y,i,i);fl=1;
  902.           if(signe(b))
  903.         {
  904.           for(k=1;(k<i)&&fl;k++)
  905.             for(l=1;(l<i)&&fl;l++)
  906.               fl= !signe(modii(coeff(y,k,l),b));
  907.         }
  908.           if(!fl)
  909.         {
  910.           k--;
  911.           for(l=1;l<=i;l++)
  912.             coeff(y,i,l)=laddii(coeff(y,i,l),coeff(y,k,l));
  913.         }
  914.         }
  915.       if(avma<lim) {tetpil=avma;y=gerepile(av,tetpil,gcopy(y));}
  916.     }
  917.       while(c||(!fl));
  918.     }
  919.   tetpil=avma;z=cgetg(n+1,17);
  920.   for(j=0,k=1;k<=n;k++) if(!signe(coeff(y,k,k))) z[++j]=zero;
  921.   for(k=1;k<=n;k++) if(signe(coeff(y,k,k))) z[++j]=labs(coeff(y,k,k));
  922.   return gerepile(av,tetpil,z);
  923. }
  924.  
  925.  
  926. GEN transroot(x,i,j)
  927.      GEN x;
  928.      int i,j;
  929. {
  930.   long n=lg(x),k;
  931.   GEN y;
  932.  
  933.   y=cgetg(n,18);
  934.   for(k=1;k<n;k++) y[k]=((k==i)||(k==j))?x[i+j-k]:x[k];
  935.   return y;
  936. }
  937.  
  938. GEN tchirnhausen(x)
  939.      GEN x;
  940. {
  941.   long av=avma,tetpil,v,n,a,b,c;
  942.   GEN u;
  943.  
  944.   if(typ(x)!=10) err(galer1);
  945.   n=lgef(x)-3;if(n<=0) err(galer1);v=varn(x);
  946.   if(v) {u=gcopy(x);setvarn(u,0);x=u;}
  947.   do
  948.     {
  949.       a=rand()&15;if(!a) a=1;b=rand()&31;if(b>=16) b-=32;
  950.       c=rand()&31;if(c>=16) c-=32;
  951.       u=caract(gmodulcp(gaddsg(c,gmul(polx[0],gaddsg(b,gmulsg(a,polx[0])))),x),v);
  952.     }
  953.   while(lgef(polgcd(u,deriv(u,v)))>=4);
  954.   tetpil=avma;return gerepile(av,tetpil,gcopy(u));
  955. }
  956.  
  957. int gpolcomp(p1,p2)
  958.      GEN p1,p2;
  959. {
  960.   int d,j;
  961.  
  962.   d=lgef(p1)-3;if((lgef(p2)-3)!=d) err(gpolcompbug1);
  963.   j=d+1;while((j>=2)&&gegal(absi(p1[j]),absi(p2[j]))) j--;
  964.   if(j==1) return 0;
  965.   return gcmp(absi(p1[j]),absi(p2[j]));
  966. }
  967.  
  968. GEN galois(x,prec)
  969.      GEN x;
  970.      long prec;
  971. {
  972.  
  973.   long av=avma,av1,i,j,k,n,f,l,l2,e,e1;
  974.   GEN x1,p1,p2,p3,p4,p5,p6,y,z;
  975.   static int ind5[20]={2,5,3,4,1,3,4,5,1,5,2,4,1,2,3,5,1,4,2,3};
  976.   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};
  977.  
  978.   if(typ(x)!=10) err(galer1);
  979.   n=lgef(x)-3;if(n<=0) err(galer1);
  980.   if(n>7) err(impl,"galois of degree higher than 7");
  981.   x=gdiv(x,content(x));
  982.   for(i=2;i<=n+2;i++) if(typ(x[i])!=1) err(galer2);
  983.   p1=(GEN)x[n+2];
  984.   if(!gcmp1(p1))
  985.     {
  986.       x1=cgetg(n+3,10);x1[1]=x[1];x1[n+2]=un;p2=gun;
  987.       for(i=n+1;i>=2;i--) {x1[i]=lmul(x[i],p2);if(i>2) p2=gmul(p1,p2);}
  988.       x=x1;
  989.     }
  990.   p1=factor(x);
  991.   if((lg(p1[1])>2)||(!gcmp1(coeff(p1,1,2)))) err(impl,"galois of reducible polynomial");
  992.   x1=gcopy(x);av1=avma;
  993.   for(;;)
  994.     {
  995.       switch(n)
  996.     {
  997.     case 1: avma=av;y=cgetg(4,17);y[1]=y[3]=un;y[2]=lneg(gun);return y;
  998.     case 2: avma=av;y=cgetg(4,17);y[1]=deux;y[3]=un;y[2]=lneg(gun);return y;
  999.     case 3: f=carreparfait(discsr(x));avma=av;y=cgetg(4,17);y[3]=un;
  1000.       if(f) {y[2]=un;y[1]=lstoi(3);return y;}
  1001.       else {y[2]=lneg(gun);y[1]=lstoi(6);return y;}
  1002.     case 4: 
  1003.       do
  1004.         {
  1005.           p1=roots(x,prec);p2=p1;
  1006.           p3=gzero;for(i=1;i<=4;i++) p3=gadd(p3,gmul(p2[i],gsqr(p2[(i&3)+1])));
  1007.           p4=gsub(polx[0],p3);p2=transroot(p1,1,2);
  1008.           p3=gzero;for(i=1;i<=4;i++) p3=gadd(p3,gmul(p2[i],gsqr(p2[(i&3)+1])));
  1009.           p4=gmul(p4,gsub(polx[0],p3));p2=transroot(p1,1,3);
  1010.           p3=gzero;for(i=1;i<=4;i++) p3=gadd(p3,gmul(p2[i],gsqr(p2[(i&3)+1])));
  1011.           p4=gmul(p4,gsub(polx[0],p3));p2=transroot(p1,1,4);
  1012.           p3=gzero;for(i=1;i<=4;i++) p3=gadd(p3,gmul(p2[i],gsqr(p2[(i&3)+1])));
  1013.           p4=gmul(p4,gsub(polx[0],p3));p2=transroot(p1,2,3);
  1014.           p3=gzero;for(i=1;i<=4;i++) p3=gadd(p3,gmul(p2[i],gsqr(p2[(i&3)+1])));
  1015.           p4=gmul(p4,gsub(polx[0],p3));p2=transroot(p1,3,4);
  1016.           p3=gzero;for(i=1;i<=4;i++) p3=gadd(p3,gmul(p2[i],gsqr(p2[(i&3)+1])));
  1017.           p4=gmul(p4,gsub(polx[0],p3));p5=grndtoi(greal(p4),&e);
  1018.           e=max(e,gexpo(gimag(p4)));
  1019.           if(e> -10) prec=(prec<<2)-2;
  1020.         }
  1021.       while(e> -10);
  1022.       p6=ggcd(p5,deriv(p5,0));
  1023.       f=(typ(p6)==10)&&(lgef(p6)>3);
  1024.       if(f) goto tchi;
  1025.       p1=factor(p5);p2=(GEN)p1[1];l=lg(p2)-1;
  1026.       switch(l)
  1027.         {
  1028.         case 1: f=carreparfait(discsr(x));avma=av;y=cgetg(4,17);y[3]=un;
  1029.           if(f) {y[2]=un;y[1]=lstoi(12);return y;}
  1030.           else {y[2]=lneg(gun);y[1]=lstoi(24);return y;}
  1031.         case 2: avma=av;y=cgetg(4,17);y[3]=un;y[2]=lneg(gun);y[1]=lstoi(8);return y;
  1032.         case 3: l2=lgef(p2[1])-3;
  1033.           if(l2==2) {avma=av;y=cgetg(4,17);y[3]=y[2]=un;y[1]=lstoi(4);return y;}
  1034.           else {avma=av;y=cgetg(4,17);y[3]=un;y[2]=lneg(gun);y[1]=lstoi(4);return y;}
  1035.         default: err(galbug1);
  1036.         }
  1037.     case 5:
  1038.       do
  1039.         {
  1040.           do
  1041.         {
  1042.           p1=roots(x,prec);z=cgetg(7,17);
  1043.           for(l=1;l<=5;l++)
  1044.             {
  1045.               p2=(l==1)?p1:transroot(p1,1,l);
  1046.               p3=gzero;k=0;for(i=1;i<=5;i++)
  1047.             {
  1048.               p5=gadd(gmul(p2[ind5[k]],p2[ind5[k+1]]),gmul(p2[ind5[k+2]],p2[ind5[k+3]]));
  1049.               p3=gadd(p3,gmul(gsqr(p2[i]),p5));k+=4;
  1050.             }
  1051.               z[l]=lrndtoi(greal(p3),&e);
  1052.               p4=(l==1)?gsub(polx[0],p3):gmul(p4,gsub(polx[0],p3));
  1053.             }
  1054.           p2=transroot(p1,2,5);
  1055.           p3=gzero;k=0;for(i=1;i<=5;i++)
  1056.             {
  1057.               p5=gadd(gmul(p2[ind5[k]],p2[ind5[k+1]]),gmul(p2[ind5[k+2]],p2[ind5[k+3]]));
  1058.               p3=gadd(p3,gmul(gsqr(p2[i]),p5));k+=4;
  1059.             }
  1060.           z[6]=lrndtoi(greal(p3),&e);
  1061.           p4=gmul(p4,gsub(polx[0],p3));
  1062.           p5=grndtoi(greal(p4),&e);
  1063.           e=max(e,gexpo(gimag(p4)));
  1064.           if(e> -10) prec=(prec<<2)-2;
  1065.         }
  1066.           while(e> -10);
  1067.           p6=ggcd(p5,deriv(p5,0));
  1068.           f=(typ(p6)==10)&&(lgef(p6)>3);
  1069.           if(f) goto tchi;
  1070.           p3=factor(p5);l=lg(p3[1])-1;
  1071.           f=carreparfait(discsr(x));
  1072.           if(l==1)
  1073.         {
  1074.           avma=av;y=cgetg(4,17);y[3]=un;
  1075.           if(f) {y[2]=un;y[1]=lstoi(60);return y;}
  1076.           else {y[2]=lneg(gun);y[1]=lstoi(120);return y;}
  1077.         }
  1078.           else
  1079.         {
  1080.           if(f) 
  1081.             {
  1082.               l=1;while((l<=6)&&(!gcmp0(poleval(p5,z[l])))) l++;
  1083.               if(l>6) err(galbug4);
  1084.               p2=(l==6)?transroot(p1,2,5):transroot(p1,1,l);
  1085.               p3=gzero;
  1086.               for(i=1;i<=5;i++)
  1087.             {
  1088.               j=(i%5)+1;
  1089.               p3=gadd(p3,gmul(gmul(p2[i],p2[j]),gsub(p2[j],p2[i])));
  1090.             }
  1091.               p5=gmul(p3,p3);p4=grndtoi(greal(p5),&e1);
  1092.               e1=max(e1,gexpo(gimag(p5)));
  1093.               if(e1<= -10)
  1094.             {
  1095.               if(gcmp0(p4)) goto tchi;
  1096.               f=carreparfait(p4);
  1097.               avma=av;y=cgetg(4,17);y[3]=y[2]=un;y[1]=lstoi(f?5:10);return y;
  1098.             }
  1099.               else prec=(prec<<2)-2;
  1100.             }
  1101.           else
  1102.             {
  1103.               avma=av;y=cgetg(4,17);y[3]=un;y[2]=lneg(gun);y[1]=lstoi(20);return y;
  1104.             }
  1105.         }
  1106.         }
  1107.       while(e1> -10);
  1108.     case 6: 
  1109.       do
  1110.         {
  1111.           do
  1112.         {
  1113.           p1=roots(x,prec);
  1114.           for(l=1;l<=6;l++)
  1115.             {
  1116.               p2=(l==1)?p1:transroot(p1,1,l);
  1117.               p3=gzero;k=0;for(i=1;i<=5;i++) for(j=i+1;j<=6;j++)
  1118.             {
  1119.               p5=gadd(gmul(p2[ind6[k]],p2[ind6[k+1]]),gmul(p2[ind6[k+2]],p2[ind6[k+3]]));
  1120.               p3=gadd(p3,gmul(gsqr(gmul(p2[i],p2[j])),p5));k+=4;
  1121.             }
  1122.               p4=(l==1)?gsub(polx[0],p3):gmul(p4,gsub(polx[0],p3));
  1123.             }
  1124.           p5=grndtoi(greal(p4),&e);
  1125.           e=max(e,gexpo(gimag(p4)));
  1126.           if(e> -10) prec=(prec<<2)-2;
  1127.         }
  1128.           while(e> -10);
  1129.           p6=ggcd(p5,deriv(p5,0));
  1130.           f=(typ(p6)==10)&&(lgef(p6)>3);
  1131.           if(f) goto tchi;
  1132.           p3=factor(p5);p2=(GEN)p3[1];l=lg(p2)-1;
  1133.           switch(l)
  1134.         {
  1135.         case 1:    p3=gadd(gmul(gmul(p1[1],p1[2]),p1[3]),gmul(gmul(p1[4],p1[5]),p1[6]));
  1136.           p4=gsub(polx[0],p3);
  1137.           for(i=1;i<=3;i++)
  1138.             for(j=4;j<=6;j++)
  1139.               {
  1140.             p2=transroot(p1,i,j);
  1141.             p3=gadd(gmul(gmul(p2[1],p2[2]),p2[3]),gmul(gmul(p2[4],p2[5]),p2[6]));
  1142.             p4=gmul(p4,gsub(polx[0],p3));
  1143.               }
  1144.           p5=grndtoi(greal(p4),&e1);
  1145.           e1=max(e1,gexpo(gimag(p4)));
  1146.           if(e1<= 10)
  1147.             {
  1148.               p6=ggcd(p5,deriv(p5,0));
  1149.               f=(typ(p6)==10)&&(lgef(p6)>3);
  1150.               if(f) goto tchi;
  1151.               p3=factor(p5);p2=(GEN)p3[1];l=lg(p2)-1;f=carreparfait(discsr(x));
  1152.               avma=av;y=cgetg(4,17);y[3]=un;
  1153.               if(l==1)
  1154.             {
  1155.               if(f) {y[2]=un;y[1]=lstoi(360);return y;}
  1156.               else {y[2]=lneg(un);y[1]=lstoi(720);return y;}
  1157.             }
  1158.               else
  1159.             {
  1160.               if(f) {y[2]=un;y[1]=lstoi(36);return y;}
  1161.               else {y[2]=lneg(un);y[1]=lstoi(72);return y;}
  1162.             }
  1163.             }
  1164.           else prec=(prec<<2)-2;
  1165.           break;
  1166.           
  1167.         case 2: l2=lgef(p2[1])-3;if(l2>3) l2=6-l2;
  1168.           switch(l2)
  1169.             {
  1170.             case 1: f=carreparfait(discsr(x));avma=av;y=cgetg(4,17);y[3]=un;
  1171.               if(f) {y[2]=un;y[1]=lstoi(60);return y;}
  1172.               else {y[2]=lneg(un);y[1]=lstoi(120);return y;}
  1173.             case 2: f=carreparfait(discsr(x));
  1174.               if(f) {avma=av;y=cgetg(4,17);y[3]=y[2]=un;y[1]=lstoi(24);return y;}
  1175.               else
  1176.             {
  1177.               p3=(lgef(p2[1])==5)?(GEN)p2[2]:(GEN)p2[1];
  1178.               f=carreparfait(discsr(p3));avma=av;y=cgetg(4,17);y[2]=lneg(gun);
  1179.               if(f) {y[1]=lstoi(24);y[3]=deux;return y;}
  1180.               else {y[1]=lstoi(48);y[3]=un;return y;}
  1181.             }
  1182.             case 3: f=carreparfait(discsr(p2[1]))||carreparfait(discsr(p2[2]));
  1183.               avma=av;y=cgetg(4,17);y[3]=un;y[2]=lneg(gun);y[1]=lstoi(f?18:36);
  1184.               return y;
  1185.             }
  1186.         case 3: for(l2=1;l2<=3;l2++) if(lgef(p2[l2])>=6) p3=(GEN)p2[l2];
  1187.           if(lgef(p3)==6)
  1188.             {
  1189.               f=carreparfait(discsr(p3));avma=av;y=cgetg(4,17);y[2]=lneg(gun);
  1190.               y[3]=un;y[1]=f?lstoi(6):lstoi(12);return y;
  1191.             }
  1192.           else
  1193.             {
  1194.               f=carreparfait(discsr(x));avma=av;y=cgetg(4,17);y[3]=un;
  1195.               if(f) {y[2]=un;y[1]=lstoi(12);return y;}
  1196.               else {y[2]=lneg(un);y[1]=lstoi(24);return y;}
  1197.             }
  1198.         case 4: avma=av;y=cgetg(4,17);y[1]=lstoi(6);y[2]=lneg(gun);
  1199.           y[3]=deux;return y;
  1200.         default: err(galbug3);
  1201.         }
  1202.         }
  1203.       while(e1> -10);
  1204.       
  1205.     case 7: 
  1206.       do
  1207.         {
  1208.           p1=roots(x,prec);p4=gun;
  1209.           for(i=1;i<=5;i++)
  1210.         for(j=i+1;j<=6;j++)
  1211.           for(k=j+1;k<=7;k++)
  1212.             p4=gmul(p4,gsub(polx[0],gadd(gadd(p1[i],p1[j]),p1[k])));
  1213.           p5=grndtoi(greal(p4),&e);e=max(e,gexpo(gimag(p4)));
  1214.           if(e> -10) prec=(prec<<2)-2;
  1215.         }
  1216.       while(e> -10);
  1217.       p6=ggcd(p5,deriv(p5,0));
  1218.       f=(typ(p6)==10)&&(lgef(p6)>3);
  1219.       if(f) goto tchi;
  1220.       p1=factor(p5);p2=(GEN)p1[1];l=lg(p2)-1;
  1221.       switch(l)
  1222.         {
  1223.         case 1: f=carreparfait(discsr(x));avma=av;y=cgetg(4,17);y[3]=un;
  1224.           if(f) {y[2]=un;y[1]=lstoi(2520);return y;}
  1225.           else {y[2]=lneg(gun);y[1]=lstoi(5040);return y;}
  1226.         case 2: f=lgef(p2[1])-3;avma=av;y=cgetg(4,17);y[3]=un;
  1227.           if((f==7)||(f==28)) {y[2]=un;y[1]=lstoi(168);return y;}
  1228.           else {y[2]=lneg(gun);y[1]=lstoi(42);return y;}
  1229.         case 3: avma=av;y=cgetg(4,17);y[3]=y[2]=un;y[1]=lstoi(21);return y;
  1230.         case 4: avma=av;y=cgetg(4,17);y[3]=un;y[2]=lneg(gun);y[1]=lstoi(14);return y;
  1231.         case 5: avma=av;y=cgetg(4,17);y[3]=y[2]=un;y[1]=lstoi(7);return y;
  1232.         default: err(galbug2);
  1233.         }
  1234.     }
  1235.     tchi:
  1236.       avma=av1;x=tchirnhausen(x1);
  1237.     }
  1238. }
  1239.  
  1240. GEN initalg(x,prec)
  1241.      GEN x;
  1242.      long prec;
  1243. {
  1244.   GEN y,p1,p2,p3,p4,p5,p6,p7,fieldd,polr,ptrace,dx,adx,polmax,s,a,dxn,adxn,sn,phimax,ind;
  1245.   long tx=typ(x),n=lgef(x)-3,i,j,av=avma,av1,v,k,r1,imax,numb,flc,tetpil;
  1246.  
  1247.   if((tx!=10)||(n<=0)) err(poltyper);
  1248.   p1=content(x);p1=gcmp1(p1) ? x: gdiv(x,content(x));
  1249.   for(k=2;k<=n+2;k++) if(typ(p1[k])!=1) err(impl,"general algebraic extension");
  1250.   if(gcmp1(p2=(GEN)p1[n+2])) x=p1;
  1251.   else
  1252.     {
  1253.       x=cgetg(n+3,10);x[1]=p1[1];x[n+2]=un;x[n+1]=p1[n+1];p3=p2;
  1254.       for(k=n;k>=2;k--)
  1255.     {
  1256.       x[k]=lmulii(p3,p1[k]);
  1257.       if(k>2) p3=mulii(p2,p3);
  1258.     }
  1259.     }
  1260.   p1=factor(x);
  1261.   if(lgef(coeff(p1,1,1))!=n+3) err(redper1);
  1262.   r1=sturm(x);p4=allbase(x,0,&fieldd);
  1263.   if(r1<n)
  1264.     {
  1265.       polr=roots(x,prec);p3=cgetg(n+1,19);
  1266.       for(i=1;i<=n;i++)
  1267.     {
  1268.       p1=cgetg(n+1,18);p3[i]=(long)p1;
  1269.       for(j=1;j<=n;j++)
  1270.         p1[j]=lsubst(p4[i],varn(p4[i]),polr[j]);
  1271.     }
  1272.       p2=greal(gmul(gconj(gtrans(p3)),p3));
  1273.     }
  1274.   else
  1275.     {
  1276.       ptrace=cgetg(n+1,17);ptrace[1]=lstoi(n);
  1277.       for(k=1;k<n;k++) 
  1278.     {
  1279.       p3=gmulsg(k,x[n-k+2]);
  1280.       for(i=1;i<k;i++) p3=gadd(p3,gmul(x[n-i+2],ptrace[k-i+1]));
  1281.       ptrace[k+1]=lneg(p3);
  1282.     }
  1283.       p2=cgetg(n+1,19);
  1284.       for(i=1;i<=n;i++)
  1285.     {
  1286.       p1=cgetg(n+1,18);p2[i]=(long)p1;
  1287.       for(j=1;j<i;j++) p1[j]=lcopy(coeff(p2,i,j));
  1288.       for(j=i;j<=n;j++)
  1289.         {
  1290.           p5=gres(gmul(p4[i],p4[j]),x);p6=gzero;
  1291.           for(k=0;k<=lgef(p5)-3;k++) p6=gadd(p6,gmul(p5[k+2],ptrace[k+1]));
  1292.           p1[j]=(long)p6;
  1293.         }
  1294.     }
  1295.     }
  1296.   p1=lllgram(p2,prec);v=varn(x);dx=discsr(x);adx=absi(dx);polmax=x;imax=0;
  1297.   if(r1<n) for(s=gzero,i=1;i<=n;i++) s=gadd(s,gnorm(polr[i]));
  1298.   else s=gsub(gsqr(x[n+1]),gmul2n(x[n],1));
  1299.   a=cgetg(n+1,18);for(i=1;i<=n;i++) a[i]=lmul(p4,p1[i]);
  1300.   for(numb=0,i=1;i<=n;i++)
  1301.     {
  1302.       av1=avma;p3=gmodulcp(a[i],x);p7=content(p3[2]);
  1303.       if(gcmp1(p7)) p3=caract(p3,v);
  1304.       else
  1305.     {
  1306.       p3=caract(gdiv(p3,p7),v);
  1307.       p3=gmul(gpuigs(p7,lgef(p3)-3),gsubst(p3,v,gdiv(polx[v],p7)));
  1308.     }
  1309.       p5=ggcd(deriv(p3,v),p3);
  1310.       if(lgef(p5)==3)
  1311.     {
  1312.       dxn=discsr(p3);adxn=absi(dxn);flc=gcmp(adxn,adx);numb++;
  1313.       if(flc<=0)
  1314.         {
  1315.           if(r1<n) for(sn=gzero,j=1;j<=n;j++) sn=gadd(sn,gnorm(poleval(a[i],polr[j])));
  1316.           else sn=gsub(gsqr(p3[n+1]),gmul2n(p3[n],1));
  1317.           if(flc<0) {dx=dxn;adx=adxn;s=sn;polmax=p3;imax=i;}
  1318.           else 
  1319.         {
  1320.           flc=gcmp(sn,s);
  1321.           if((flc<0)||((!flc)&&(gpolcomp(p3,polmax)<0)))
  1322.             {dx=dxn;adx=adxn;s=sn;polmax=p3;imax=i;}
  1323.         }
  1324.         }
  1325.     }
  1326.     }
  1327.   if(!numb) err(polreder1);
  1328.   phimax=imax?(GEN)a[imax]:polx[v];
  1329.   j=n+1;
  1330.   while((j>=2)&&(!signe(polmax[j]))) j-=2;
  1331.   if((j>=2)&&(signe(polmax[j])>0))
  1332.     {
  1333.       for(;j>=2;j-=2) setsigne(polmax[j],-signe(polmax[j]));
  1334.       phimax=gneg(phimax);
  1335.     }
  1336.   p2=lift(gsubst(p4,v,polymodrecip(gmodulcp(phimax,x))));
  1337.   p3=cgetg(n+1,19);
  1338.   for(j=1;j<=n;j++)
  1339.     {
  1340.       p4=cgetg(n+1,18);p3[j]=(long)p4;
  1341.       for(i=1;i<=n;i++) p4[i]=(long)truecoeff(p2[j],i-1);
  1342.     }
  1343.   p4=denom(p3);
  1344.   p2=gdiv(hnf(gmul(p4,p3)),p4);p3=cgetg(n+1,17);
  1345.   for(j=1;j<=n;j++) 
  1346.     {
  1347.       p1=gzero;for(i=n;i;i--) p1=gadd(coeff(p2,i,j),gmul(p1,polx[v]));
  1348.       p3[j]=(long)p1;
  1349.     }
  1350.   if(!carrecomplet(divii(dx,fieldd),&ind)) err(initalgbug1);
  1351.   tetpil=avma;y=cgetg(8,17);y[1]=lcopy(polmax);p1=cgetg(3,17);
  1352.   p1[1]=lstoi(r1);p1[2]=lstoi((n-r1)>>1);y[2]=(long)p1;
  1353.   y[3]=lcopy(fieldd);y[4]=lcopy(ind);y[5]=lcopy(s);
  1354.   if(r1<n) 
  1355.     {
  1356.       p1=cgetg(n+1,18);
  1357.       for(i=1;i<=n;i++) p1[i]=(long)poleval(phimax,polr[i]);
  1358.       y[6]=(long)p1;
  1359.     }
  1360.   else y[6]=lreal(roots(polmax,prec));
  1361.   y[7]=lcopy(p3);
  1362.   return gerepile(av,tetpil,y);
  1363. }
  1364.  
  1365.  
  1366. GEN galoisconj2(x,prec)
  1367.      GEN x;
  1368.      long prec;
  1369. {
  1370.   long av=avma,tetpil,i,j,n,v;
  1371.   GEN y,w,polr,p1,p2,b,di;
  1372.  
  1373.   if(typ(x)!=10) err(galer1);
  1374.   n=lgef(x)-3;if(n<=0) err(galer1);
  1375.   v=varn(x);p1=factor(x);
  1376.   if((lg(p1[1])>2)||(!gcmp1(coeff(p1,1,2)))) err(galcer1);
  1377.   polr=roots(x,prec);p1=(GEN)polr[1];b=allbase(x,0,&di);
  1378.   w=cgetg(n+1,17);w[1]=un;for(i=2;i<=n;i++) w[i]=(long)gsubst(b[i],v,p1);
  1379.   y=cgetg(n+1,17);y[1]=(long)polx[v];
  1380.   for(i=2;i<=n;i++)
  1381.     {
  1382.       p1=lindep(concat(w,polr[i]),prec);
  1383.       if(gcmp0(p1[n+1])) y[i]=zero;
  1384.       else
  1385.     {
  1386.       p2=gzero;
  1387.       for(j=1;j<=n;j++) p2=gadd(p2,gmul(p1[j],b[j]));
  1388.       p2=gdiv(p2,gneg(p1[n+1]));
  1389.       if(gcmp0(gmod(gsubst(x,v,p2),x))) y[i]=(long)p2;else y[i]=zero;
  1390.     }
  1391.     }
  1392.   tetpil=avma;return gerepile(av,tetpil,gcopy(y));
  1393. }
  1394.  
  1395. GEN galoisconj(x,prec)
  1396.      GEN x;
  1397.      long prec;
  1398. {
  1399.   long av=avma,tetpil,i,j,n,v;
  1400.   GEN y,w,polr,p1,p2;
  1401.  
  1402.   if(typ(x)!=10) err(galer1);
  1403.   n=lgef(x)-3;if(n<=0) err(galer1);
  1404.   v=varn(x);p1=factor(x);
  1405.   if((lg(p1[1])>2)||(!gcmp1(coeff(p1,1,2)))) err(galcer1);
  1406.   polr=roots(x,prec);p1=(GEN)polr[1];
  1407.   w=cgetg(n+1,17);w[1]=un;for(i=2;i<=n;i++) w[i]=lmul(p1,w[i-1]);
  1408.   y=cgetg(n+1,17);y[1]=(long)polx[v];
  1409.   for(i=2;i<=n;i++)
  1410.     {
  1411.       p1=lindep(concat(w,polr[i]),prec);
  1412.       if(gcmp0(p1[n+1])) y[i]=zero;
  1413.       else
  1414.     {
  1415.       p2=gzero;
  1416.       for(j=n;j;j--) p2=gadd(p1[j],gmul(p2,polx[v]));
  1417.       p2=gdiv(p2,gneg(p1[n+1]));
  1418.       if(gcmp0(gmod(gsubst(x,v,p2),x))) y[i]=(long)p2;else y[i]=zero;
  1419.     }
  1420.     }
  1421.   tetpil=avma;return gerepile(av,tetpil,gcopy(y));
  1422. }
  1423.  
  1424.  
  1425.  
  1426.  
  1427.  
  1428.  
  1429.  
  1430.  
  1431.  
  1432.  
  1433.