home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Professional / OS2PRO194.ISO / os2 / progs / pari / pari_137 / src / bibli1.c < prev    next >
Text File  |  1992-05-20  |  68KB  |  2,648 lines
  1. /********************************************************************/
  2. /********************************************************************/
  3. /**                                                                **/
  4. /**                  BIBLIOTHEQUE  MATHEMATIQUE                    **/
  5. /**                     (premiere partie)                          **/
  6. /**                                                                **/
  7. /**                     Copyright Babe Cool                        **/
  8. /**                                                                **/
  9. /********************************************************************/
  10. /********************************************************************/
  11.  
  12. # include "genpari.h"
  13.  
  14.  
  15. /********************************************************************/
  16. /********************************************************************/
  17. /**                                                                **/
  18. /**                 DEVELOPPEMENTS  LIMITES                        **/
  19. /**                                                                **/
  20. /********************************************************************/
  21. /********************************************************************/
  22.  
  23. GEN     tayl(x,v,precdl)
  24.      
  25.      GEN        x;
  26.      long       v, precdl;
  27.      
  28. {
  29.   long tetpil,i,vx = gvar9(x), av=avma;
  30.   GEN p1, y;
  31.   
  32.   if(v <= vx)
  33.   {
  34.     p1 = cgetg(3,11);
  35.     p1[2] = zero; p1[1] = 0x8000 + precdl; setvarn(p1,v);
  36.     tetpil = avma; y = gerepile(av, tetpil, gadd(p1,x));
  37.   }
  38.   else
  39.   {
  40.     p1=cgetg(v+2,17);
  41.     for(i=0;i<vx;i++) p1[i+1]=lpolx[i];p1[vx+1]=lpolx[v];
  42.     for(i=vx+1;i<v;i++) p1[i+1]=lpolx[i];p1[v+1]=lpolx[vx];
  43.     y = tayl(changevar(x, p1), vx, precdl); tetpil=avma;
  44.     y = gerepile(av, tetpil, changevar(y, p1));
  45.   }
  46.   return y;
  47. }
  48.  
  49.  
  50. GEN     ggrando(x, n)
  51.      
  52.      GEN        x;
  53.      long       n;
  54.      
  55. {
  56.   long m, v, tx=typ(x);
  57.   GEN y;
  58.   if (gcmp0(x)) err(grandoer1);
  59.   if (tx>9)
  60.   {
  61.     y=cgetg(3,11);m=gval(x,v = gvar(x));
  62.     y[1]=0x8000+n*m;y[2]=zero;setvarn(y,v);
  63.   }
  64.   else
  65.   {
  66.     if(tx!=1) err(grandoer1);
  67.     y=cgetg(5,7);y[2]=lclone(x);
  68.     y[3]=un;y[4]=zero;
  69.     setvalp(y,n);setprecp(y,0);
  70.   }
  71.   return y;
  72. }
  73.  
  74. /********************************************************************/
  75. /********************************************************************/
  76. /**                                                                **/
  77. /**                POLYNOMES  DE  TCHEBICHEFF                      **/
  78. /**                                                                **/
  79. /********************************************************************/
  80. /********************************************************************/
  81.  
  82. /*      T0=1; T1=X ;  T(n)=2*X*T(n-1)-T(n-2)    */
  83.  
  84. GEN  tchebi(n)
  85.      
  86.      long  n;
  87.      
  88. {
  89.   long  av1,av2,av3,av4,av5,decal,m;
  90.   GEN  p0,p1,p2,px,q;
  91.   
  92.   if (n==0) return polun[0];
  93.   else
  94.     {
  95.       av1=avma;
  96.       px=cgetg(4,10);
  97.       px[1]=0x1000004;px[2]=zero;px[3]=un;
  98.       if (n==1) return px;
  99.       else
  100.         {
  101.           av2=avma;
  102.           p0=gcopy(polun[0]);
  103.           av3=avma;
  104.           p1=gcopy(px);
  105.           for (m=1;m<n-1;m++)
  106.             {
  107.               q=gmul(px,gmulsg(2,p1));
  108.               av4=avma;
  109.               p2=gsub(q,p0);av5=avma;
  110.               decal=lpile(av2,av3,0);
  111.               p0=adecaler(p1,av3,av5)?p1+(decal>>2):p1;
  112.           p1=adecaler(p2,av3,av5)?p2+(decal>>2):p2;
  113.               av3=av4+decal;
  114.             }
  115.           q=gmul(px,gmulsg(2,p1));av2=avma;
  116.           return gerepile(av1,av2,gsub(q,p0));
  117.         }
  118.     }
  119. }
  120.  
  121. /********************************************************************/
  122. /********************************************************************/
  123. /**                                                                **/
  124. /**                  POLYNOMES  DE  LEGENDRE                       **/
  125. /**                                                                **/
  126. /********************************************************************/
  127. /********************************************************************/
  128.  
  129. /* L0=1; L1=X ;(n+1)*L(n)=(2*n+1)*X*L(n)-n*L(n-1) */
  130.  
  131. GEN  legendre(n)
  132.      
  133.      long  n;
  134.      
  135. {
  136.   long  av1,av2,av3,av4,av5,av6,m,decal;
  137.   GEN  p0,p1,p2,px,q2;
  138.   
  139.   if (n==0) return polun[0];
  140.   else
  141.   {
  142.     av1=avma;
  143.     px=cgetg(4,10);
  144.     px[1]=0x1000004;px[2]=zero;px[3]=un;
  145.     if (n==1) return px;
  146.     else
  147.     {
  148.       av2=avma;
  149.       p0=polun[0];
  150.       av3=avma;
  151.       p1=gmulsg(2,px);
  152.       
  153.       for (m=1;m<n;m++)
  154.       {
  155.         av4=avma;
  156.         p2=gsub(gmul(px,gmulsg(4*m+2,p1)),gmulsg(4*m,p0));
  157.         av5=avma;
  158.         p2=gerepile(av4,av5,gdivgs(p2,m+1));
  159.         av6=avma;decal=lpile(av2,av3,0);
  160.         p0=adecaler(p1,av3,av6)?p1+(decal>>2):p1;
  161.     p1=adecaler(p2,av3,av6)?p2+(decal>>2):p2;
  162.         av3=av4+decal;
  163.       }
  164.       q2=mpshift(un,n);av2=avma;
  165.       return gerepile(av1,av2,gdiv(p1,q2));
  166.     }
  167.   }
  168. }
  169.  
  170. GEN cyclo(n)
  171.      long n;
  172. {
  173.   long av=avma,tetpil,d,q,m;
  174.   GEN p1,yn,yd;
  175.  
  176.   if(n<=0) err(arither2);
  177.   d=1;yn=gun;yd=gun;
  178.   while(d*d<=n)
  179.     {
  180.       if(!(n%d))
  181.     {
  182.       q=n/d;m=mu(stoi(q));
  183.       if(m)
  184.         {
  185.           p1=gsub(gpuigs(polx[0],d),gun);
  186.           if(m>0) yn=gmul(yn,p1);else yd=gmul(yd,p1);
  187.         }
  188.       if(q!=d)
  189.         {
  190.           m=mu(stoi(d));
  191.           if(m)
  192.         {
  193.           p1=gsub(gpuigs(polx[0],q),gun);
  194.           if(m>0) yn=gmul(yn,p1);else yd=gmul(yd,p1);
  195.         }
  196.         }
  197.     }
  198.       d++;
  199.     }
  200.   tetpil=avma;return gerepile(av,tetpil,gdiv(yn,yd));
  201. }
  202.           
  203.  
  204. /********************************************************************/
  205. /********************************************************************/
  206. /**                                                                **/
  207. /**              MATRICES DE HILBERT,PASCAL ...                    **/
  208. /**                                                                **/
  209. /********************************************************************/
  210. /********************************************************************/
  211.  
  212. GEN     hilb(n)
  213.      /* matrice de hilbert d'ordre n */
  214.      long    n;
  215.      
  216.      
  217. {
  218.   long i,j,sho[3];GEN p,a;
  219.   p=cgetg(n+1,19);
  220.   sho[0] = 0x1010003;
  221.   sho[1] = 0x1000003;
  222.   for (j=1;j<=n;j++)
  223.   {
  224.     p[j]=lgetg(n+1,18);
  225.     for (i=1;i<=n;i++)
  226.     {
  227.       a=cgetg(3,4);coeff(p,i,j)=(long)a;
  228.       a[1]=un;sho[2]=i+j-1;a[2]=lcopy(sho);
  229.     }
  230.   }
  231.   return p;
  232. }
  233.  
  234. GEN     pasc(n)
  235.      /* matrice triangle de PASCAL d'ordre n */
  236.      long    n;
  237.      
  238. {
  239.   long i,j;GEN p;
  240.   
  241.   p=cgetg(n+2,19);
  242.   for (j=1;j<=n+1;j++)  p[j]=lgetg(n+2,18);
  243.   for (i=1;i<=n+1;i++)
  244.   {
  245.     coeff(p,i,1)=un;coeff(p,i,i)=un;
  246.     for (j=2;j<i;j++)
  247.       coeff(p,i,j)=laddii(coeff(p,i-1,j),coeff(p,i-1,j-1));
  248.     for (j=i+1;j<=n+1;j++)
  249.       coeff(p,i,j)=zero;
  250.   }
  251.   return p;
  252. }
  253.  
  254. /********************************************************************/
  255. /********************************************************************/
  256. /**                                                                **/
  257. /**              TRANSFORMEE DE LAPLACE D UNE SERIE                **/
  258. /**                                                                **/
  259. /********************************************************************/
  260. /********************************************************************/
  261.  
  262. GEN    laplace(x)
  263.      
  264.      GEN    x;
  265.      
  266. {
  267.   
  268.   long i,l,dec,ec,av1,av2,av3,av4;
  269.   GEN  y,p1,p2;
  270.   
  271.   if(typ(x)!=11) err(laper1);
  272.   if(gcmp0(x)) y=gcopy(x);
  273.   else
  274.   {
  275.     ec=valp(x);if (ec<0) err(laper2);
  276.     l=lg(x);y=cgetg(l,11);av1=avma;
  277.     p1=mpfact(ec);y[1]=x[1];
  278.     for(i=2;i<l;i++)
  279.     {
  280.       av2=avma;p2=gmul(p1,x[i]);
  281.       av3=avma;ec++;p1=gmulsg(ec,p1);
  282.       av4=avma;dec=lpile(av1,av2,0)>>2;
  283.       y[i]=adecaler(p2,av2,av4)?(long)(p2+dec):(long)p2;
  284.       if(adecaler(p1,av2,av4)) p1+=dec;
  285.       av1=av3+(dec<<2);
  286.     }
  287.     avma=av1;
  288.   }
  289.   return y;
  290. }
  291.  
  292. /********************************************************************/
  293. /********************************************************************/
  294. /**                                                                **/
  295. /**            PRODUIT DE CONVOLUTION DE DEUX SERIES               **/
  296. /**                                                                **/
  297. /********************************************************************/
  298. /********************************************************************/
  299.  
  300. GEN    convol(x,y)
  301.      
  302.      GEN    x,y;
  303.      
  304. {
  305.   long lx,ly,l,ex,ey,l1,i,j,v,f;
  306.   GEN  z;
  307.   
  308.   if((typ(x)!=11) || (typ(y)!=11)) err(convol1);
  309.   if(gcmp0(x) || gcmp0(y)) err(convol2);
  310.   lx=lg(x);ly=lg(y);ex=valp(x);ey=valp(y);
  311.   v=ex;if(ey>v) v=ey;
  312.   l=ex+lx;l1=ey+ly;if(l1<l) l=l1;
  313.   l=l-v;if(l<3) err(convol3);f=1;
  314.   for(i=v+2;(i<(l+v))&&f;i++) f=(gcmp0(x[i-ex]) || gcmp0(y[i-ey]));
  315.   if(i==(l+v)) 
  316.   {
  317.     z=cgetg(3,11); z[1]=0x7ffe +l+v; z[2]=zero;
  318.   }
  319.   else
  320.   {
  321.     z=cgetg(l-i+3+v,11); z[1]=0x1007ffd+i;
  322.     for(j=i-1;j<l+v;j++) z[j-i+3]=lmul(x[j-ex],y[j-ey]);
  323.   }
  324.   return z;
  325. }
  326.  
  327. /********************************************************************/
  328. /********************************************************************/
  329. /**                                                                **/
  330. /**         Conversion serie----->polynome ou fr. rat.             **/
  331. /**                                                                **/
  332. /**                            et                                  **/
  333. /**                                                                **/
  334. /**               p-adique-->entier ou rationnel                   **/
  335. /**                                                                **/
  336. /********************************************************************/
  337. /********************************************************************/
  338.  
  339. GEN     gconvsp(x)
  340.      
  341.      GEN x;
  342.      
  343. {
  344.   long  e ,av ,tetpil ,i ,v;
  345.   GEN   p1 ,y;
  346.   
  347.   if(typ(x)!=11) err(csper1);
  348.   v=varn(x);e=valp(x);av=avma;y=gcopy(x);settyp(y ,10);
  349.   if (gcmp0(x)) {y[1]=2;setvarn(y,v);}
  350.   else
  351.     {
  352.       for(i=lg(x)-1;(i>1)&&gcmp0(y[i]);--i);
  353.       y[1]=0x1000001+i;setvarn(y,v);
  354.       p1=gpuigs(polx[v],e);tetpil=avma;
  355.       y=gerepile(av ,tetpil ,gmul(p1 ,y));
  356.     }
  357.   return y;
  358. }
  359.  
  360. GEN gconvpe(x)
  361.      GEN x;
  362.  
  363. {
  364.   long av,tetpil;
  365.   GEN p1;
  366.  
  367.   if(typ(x)!=7) err(csper1);
  368.   if(!signe(x[4])) return gzero;
  369.   av=avma;p1=gpuigs(x[2],valp(x));
  370.   tetpil=avma;return gerepile(av,tetpil,gmul(p1,x[4]));
  371. }
  372.  
  373. /********************************************************************/
  374. /********************************************************************/
  375. /**                                                                **/
  376. /**             CHANGEMENT DE PRECISION D'UN GEN                   **/
  377. /**                                                                **/
  378. /**                       OU DU DEFAUT                             **/
  379. /**                                                                **/
  380. /********************************************************************/
  381. /********************************************************************/
  382.  
  383. GEN     gprec(x,l)
  384.      
  385.      GEN        x;
  386.      long       l;
  387.      
  388. {
  389.   long  tx=typ(x),lx=lg(x),i,pr;
  390.   GEN   y;
  391.   
  392.   if(l<=0) err(precer1);
  393.   switch(tx)
  394.     {
  395.     case 2 : pr=l*K1+3;y=cgetr(pr);affrr(x,y);break;
  396.     case 7 : y=cgetg(lx,tx);y[1]=x[1];
  397.       setprecp(y,l);y[2]=x[2];
  398.       y[3]=lpuigs(x[2],l);
  399.       y[4]=lmodii(x[4],y[3]);break;
  400.     case 11: 
  401.       if(gcmp0(x))
  402.         {
  403.           y=cgetg(3,11);y[1]=0x8000+l;setvarn(y,varn(x));
  404.         }
  405.       else
  406.         {
  407.           y=cgetg(l+2,11); y[1]=x[1];
  408.           if(l+1<lx)
  409.             {
  410.               for(i=2;i<l+2;i++)
  411.                 y[i]=lcopy(x[i]);
  412.             }
  413.           else
  414.             {
  415.               for(i=2;i<lx;i++)
  416.                 y[i]=lcopy(x[i]);
  417.               for(i=lx;i<l+2;i++)
  418.                 y[i]=zero;
  419.             }
  420.         }
  421.       break;
  422.     case 6 :
  423.     case 13:
  424.     case 14:
  425.     case 17:
  426.     case 18:
  427.     case 19: y=cgetg(lx,tx);
  428.       for(i=1;i<lx;i++)
  429.         y[i]=lprec(x[i],l);
  430.       break;
  431.     default: y=gcopy(x);
  432.     }
  433.   return y;
  434. }
  435.  
  436. long setprecr(n)
  437.      long n;
  438. {
  439.   long m=(prec-3)/K1;
  440.   if(n>0) {glbfmt[2] = n;prec = n * K1 + 3;}
  441.   return m;
  442. }
  443.  
  444. long setserieslength(n)
  445.      long n;
  446. {
  447.   long m=precdl;
  448.   if(n>0) precdl=n;
  449.   return m;
  450. }
  451.  
  452. /********************************************************************/
  453. /********************************************************************/
  454. /**                                                                **/
  455. /**                     POLYNOME RECIPROQUE       .                **/
  456. /**                                                                **/
  457. /********************************************************************/
  458. /********************************************************************/
  459.  
  460. GEN     polrecip(x)
  461.      
  462.      GEN        x;
  463.      
  464. {
  465.   GEN   y;
  466.   long  lx=lgef(x),i;
  467.   
  468.   if(typ(x)!=10) err(polrecer1);
  469.   y=cgetg(lx,10);y[1]=x[1];
  470.   for(i=2;i<lx;i++)
  471.     y[i]=lcopy(x[lx+1-i]);
  472.   normalizepol(&y);return y;
  473. }
  474.  
  475. /********************************************************************/
  476. /********************************************************************/
  477. /**                                                                **/
  478. /**                     COEFFICIENTS BINOMIAUX    .                **/
  479. /**                                                                **/
  480. /********************************************************************/
  481. /********************************************************************/
  482.  
  483. GEN     binome(x ,k)
  484.      
  485.      GEN        x;
  486.      long       k;
  487.      
  488. {
  489.   GEN   y,p1;
  490.   long  av,tetpil,i;
  491.   
  492.   if (k<0) return gzero;
  493.   if (!k) return gun;
  494.   av=avma;y=gcopy(x);if(k==1) return y;
  495.   for(i=k;i>=2;i--)
  496.     {
  497.       p1=gdivgs(gaddgs(x,i-1-k),i);tetpil=avma;
  498.       y=gmul(y,p1);
  499.     }
  500.   return gerepile(av,tetpil,y);
  501. }
  502.  
  503. /********************************************************************/
  504. /********************************************************************/
  505. /**                                                                **/
  506. /**                         ALGORITHME LLL       .                 **/
  507. /**                                                                **/
  508. /********************************************************************/
  509. /********************************************************************/
  510.  
  511. #define cosp(i,j)   (((GEN)(p3[i]))[j])
  512.  
  513. GEN    gscal(x,y)
  514.      
  515.      GEN  x,y;
  516.      
  517.      /* produit scalaire de x par y sans verification de compatibilite */
  518.      
  519. {
  520.   GEN  z,p;
  521.   long i,av,tetpil,lx=lg(x);
  522.   
  523.   z=gzero;av=avma;
  524.   for (i=1;i<lx;i++)
  525.     {
  526.       p=gmul(x[i],y[i]);
  527.       tetpil=avma;
  528.       z=gadd(z,p);
  529.     }
  530.   return gerepile(av,tetpil,z);
  531. }
  532.  
  533. GEN    lll1(x,prec)
  534.      GEN   x;
  535.      long  prec;
  536.  
  537.     /*     x est la matrice d'une base bi ;
  538.        retourne la matrice u ( entiere unimodulaire )
  539.            d'une base LLL-reduite ci en fonction des bi.
  540.       ( la base reduite est  c=b*u )                         */
  541. {
  542.   long lx=lg(x), tx=typ(x),i,j,av,av1;
  543.   GEN g;
  544.  
  545.   if(tx!=19) err(lller1);
  546.   av=avma;
  547.   g=cgetg(lx,19);
  548.   for(j=1;j<lx;j++) g[j]=lgetg(lx,18);
  549.   for(i=1;i<lx;i++)
  550.     for(j=1;j<=i;j++) coeff(g,i,j)=coeff(g,j,i)=(long)gscal(x[i],x[j]);
  551.   av1=avma;
  552.   return gerepile(av,av1,lllgram1(g,prec));
  553. }
  554.  
  555.  
  556. GEN    lllgram(x,prec)
  557.      GEN   x;
  558.      long  prec;     /* x est ici la matrice de GRAM des bi.  */
  559.      
  560. {
  561.   
  562.   GEN  mu,u,B,BB,BK,p,q,r,cst,unreel,sv;
  563.   long av0,av,tetpil,lim,l,i,j,k,lx=lg(x),tx=typ(x),n,temp,mu1,mu2,e;
  564.  
  565.   if(tx!=19) err(lllger1);
  566.   if(lg(x[1])!=lx) err(lllger2);n=lx-1;
  567.   if(n<=1) return idmat(n);
  568.   av0=avma;
  569.   cst=gdivgs(stoi(99),100); /* LLL proposent 0.75 */
  570.   lim=(avma+bot)>>1;
  571.   if (prec)
  572.     {
  573.       affsr(1,unreel=cgetr(prec+1));
  574.       x=gmul(x,unreel);
  575.       cst=gmul(cst,unreel);
  576.     }
  577.   av=avma;
  578.   mu=sqred3(x);B=cgetg(lx,18);
  579. /*  output(mu); */
  580.   for(i=1,l=0;i<=n;i++)
  581.     {
  582.       if(gsigne(B[i]=coeff(mu,i,i))>0) l++;
  583.       coeff(mu,i,i)=un;
  584.     }
  585.   if(l<n) err(lllger3);
  586.  
  587.   u=idmat(n);
  588.   k=2;
  589.   do
  590.     {
  591. /*      printf(" %d",k);fflush(stdout); */
  592.       if(!gcmp0(r=grndtoi(coeff(mu,k,k-1),&e)))
  593.         {
  594. /*      printf("k, l, e, r: %d, %d, %d, ",k,k-1,e);output(r); */
  595.           u[k]=lsub(u[k],gmul(r,u[k-1]));
  596.           for(j=1;j<k-1;j++)
  597.             coeff(mu,k,j)=lsub(coeff(mu,k,j),gmul(r,coeff(mu,k-1,j)));
  598.           mu1=coeff(mu,k,k-1)=lsub(coeff(mu,k,k-1),r);
  599.         }
  600.       else mu1=coeff(mu,k,k-1);
  601.       q=gmul(B[k-1],gsub(cst,mu2=lsqr(mu1)));
  602.       if(gcmp(q,B[k])>0)
  603.         {
  604.           BB=gadd(B[k],gmul(B[k-1],mu2));
  605.           coeff(mu,k,k-1)=ldiv(gmul(mu1,B[k-1]),BB);
  606.           B[k]=lmul(B[k-1],BK=gdiv(B[k],BB));
  607.           B[k-1]=(long)BB;
  608.           temp=u[k];u[k]=u[k-1];u[k-1]=temp;
  609.           for(j=1;j<=k-2;j++)
  610.             {
  611.               temp=coeff(mu,k,j);coeff(mu,k,j)=coeff(mu,k-1,j);
  612.               coeff(mu,k-1,j)=temp;
  613.             }
  614.           for(i=k+1;i<=n;i++)
  615.             {
  616.               p=(GEN)coeff(mu,i,k);
  617.           coeff(mu,i,k)=lsub(coeff(mu,i,k-1),gmul(mu1,p));
  618.               coeff(mu,i,k-1)=ladd(gmul(BK,p),gmul(coeff(mu,k,k-1),coeff(mu,i,k-1)));
  619.             }
  620.           if(k>2) k--;
  621.         }
  622.       else
  623.         {
  624.           for(l=k-2;l;l--)
  625.             {
  626.               if(!gcmp0(r=grndtoi(coeff(mu,k,l),&e)))
  627.                 {
  628. /*          printf("k, l, e, r: %d, %d, %d, ",k,l,e);output(r); */
  629.                   u[k]=lsub(u[k],gmul(r,u[l]));
  630.                   for(j=1;j<l;j++)
  631.                     coeff(mu,k,j)=lsub(coeff(mu,k,j),gmul(r,coeff(mu,l,j)));
  632.                   coeff(mu,k,l)=lsub(coeff(mu,k,l),r);
  633.                 }
  634.             }
  635.           k++;
  636.         }
  637.       if(avma<lim)
  638.     {
  639.       tetpil=avma;
  640.       sv=cgetg(4,17);
  641.       sv[1]=lcopy(B);sv[2]=lcopy(u);sv[3]=lcopy(mu);
  642.       sv=gerepile(av,tetpil,sv);
  643.       B=(GEN)sv[1];u=(GEN)sv[2];mu=(GEN)sv[3];
  644.     }
  645.     }
  646.   while(k<=n);
  647. /*  printf("\n"); */
  648.   tetpil=avma;return gerepile(av,tetpil,gcopy(u));
  649. }
  650.  
  651. GEN    lll(x,prec)
  652.      GEN   x;
  653.      long  prec;
  654.  
  655.     /*     x est la matrice d'une base bi ;
  656.        retourne la matrice u ( entiere unimodulaire )
  657.            d'une base LLL-reduite ci en fonction des bi.
  658.       ( la base reduite est  c=b*u )                         */
  659. {
  660.   long lx=lg(x), tx=typ(x),i,j,av,av1;
  661.   GEN g;
  662.  
  663.   if(tx!=19) err(lller1);
  664.   av=avma;
  665.   g=cgetg(lx,19);
  666.   for(j=1;j<lx;j++) g[j]=lgetg(lx,18);
  667.   for(i=1;i<lx;i++)
  668.     for(j=1;j<=i;j++) coeff(g,i,j)=coeff(g,j,i)=(long)gscal(x[i],x[j]);
  669.   av1=avma;
  670.   return gerepile(av,av1,lllgram(g,prec));
  671. }
  672.  
  673. GEN    lllint(x)
  674.      GEN   x;
  675.  
  676. {
  677.   long lx=lg(x), tx=typ(x),i,j,av,av1;
  678.   GEN g;
  679.  
  680.   if(tx!=19) err(lller1);
  681.   av=avma;
  682.   g=cgetg(lx,19);
  683.   for(j=1;j<lx;j++) g[j]=lgetg(lx,18);
  684.   for(i=1;i<lx;i++)
  685.     for(j=1;j<=i;j++) coeff(g,i,j)=coeff(g,j,i)=(long)gscal(x[i],x[j]);
  686.   av1=avma;return gerepile(av,av1,lllgramall(g,2));
  687. }
  688.  
  689. GEN lllkerim(x)
  690.      GEN x;
  691. {
  692.   long lx=lg(x), tx=typ(x),i,j,av,av1;
  693.   GEN g;
  694.  
  695.   if(tx!=19) err(lller1);
  696.   av=avma;
  697.   g=cgetg(lx,19);
  698.   for(j=1;j<lx;j++) g[j]=lgetg(lx,18);
  699.   for(i=1;i<lx;i++)
  700.     for(j=1;j<=i;j++) coeff(g,i,j)=coeff(g,j,i)=(long)gscal(x[i],x[j]);
  701.   av1=avma;return gerepile(av,av1,lllgramall(g,0));
  702. }
  703.  
  704. GEN    lllrat(x)
  705.      GEN   x;
  706.      
  707. {
  708.   return lll(x,0);
  709. }
  710.  
  711. GEN    lllgram1(x,prec)
  712.      GEN   x;
  713.      long  prec;     /* x est ici la matrice de GRAM des bi.  */
  714.      
  715. {
  716.   
  717.   GEN  mu,u,B,BB,BK,p,q,r,cst,unreel,sv;
  718.   long av0,av,tetpil,lim,l,i,j,k,lx=lg(x),tx=typ(x),n,temp,mu1,mu2,e;
  719.  
  720.   if(tx!=19) err(lllger1);
  721.   if(lg(x[1])!=lx) err(lllger2);n=lx-1;
  722.   if(n<=1) return idmat(n);
  723.   av0=avma;
  724.   cst=gdivgs(stoi(99),100); /* LLL proposent 0.75 */
  725.   lim=(avma+bot)>>1;
  726.   if (prec)
  727.     {
  728.       affsr(1,unreel=cgetr(prec+1));
  729.       x=gmul(x,unreel);
  730.       cst=gmul(cst,unreel);
  731.     }
  732.   av=avma;
  733.   mu=gtrans(sqred(x)); B=cgetg(lx,18);
  734.   for(i=1,l=0;i<=n;i++)
  735.     {
  736.       if(gsigne(B[i]=coeff(mu,i,i))>0) l++;
  737.       coeff(mu,i,i)=un;
  738.     }
  739.   if(l<n) err(lllger3);
  740.  
  741.   u=idmat(n);
  742.   k=2;
  743.   do
  744.     {
  745.       if(!gcmp0(r=grndtoi(coeff(mu,k,k-1),&e)))
  746.         {
  747.           u[k]=lsub(u[k],gmul(r,u[k-1]));
  748.           for(j=1;j<k-1;j++)
  749.             coeff(mu,k,j)=lsub(coeff(mu,k,j),gmul(r,coeff(mu,k-1,j)));
  750.           mu1=coeff(mu,k,k-1)=lsub(coeff(mu,k,k-1),r);
  751.         }
  752.       else mu1=coeff(mu,k,k-1);
  753.       q=gmul(B[k-1],gsub(cst,mu2=lsqr(mu1)));
  754.       if(gcmp(q,B[k])>0)
  755.         {
  756.           BB=gadd(B[k],gmul(B[k-1],mu2));
  757.           coeff(mu,k,k-1)=ldiv(gmul(mu1,B[k-1]),BB);
  758.           B[k]=lmul(B[k-1],BK=gdiv(B[k],BB));
  759.           B[k-1]=(long)BB;
  760.           temp=u[k];u[k]=u[k-1];u[k-1]=temp;
  761.           for(j=1;j<=k-2;j++)
  762.             {
  763.               temp=coeff(mu,k,j);coeff(mu,k,j)=coeff(mu,k-1,j);
  764.               coeff(mu,k-1,j)=temp;
  765.             }
  766.           for(i=k+1;i<=n;i++)
  767.             {
  768.               p=(GEN)coeff(mu,i,k);
  769.           coeff(mu,i,k)=lsub(coeff(mu,i,k-1),gmul(mu1,p));
  770.               coeff(mu,i,k-1)=ladd(gmul(BK,p),gmul(coeff(mu,k,k-1),coeff(mu,i,k-1)));
  771.             }
  772.           if(k>2) k--;
  773.         }
  774.       else
  775.         {
  776.           for(l=k-2;l;l--)
  777.             {
  778.               if(!gcmp0(r=grndtoi(coeff(mu,k,l),&e)))
  779.                 {
  780.                   u[k]=lsub(u[k],gmul(r,u[l]));
  781.                   for(j=1;j<l;j++)
  782.                     coeff(mu,k,j)=lsub(coeff(mu,k,j),gmul(r,coeff(mu,l,j)));
  783.                   coeff(mu,k,l)=lsub(coeff(mu,k,l),r);
  784.                 }
  785.             }
  786.           k++;
  787.         }
  788.       if(avma<lim)
  789.     {
  790.       tetpil=avma;
  791.       sv=cgetg(4,17);
  792.       sv[1]=lcopy(B);sv[2]=lcopy(u);sv[3]=lcopy(mu);
  793.       sv=gerepile(av,tetpil,sv);
  794.       B=(GEN)sv[1];u=(GEN)sv[2];mu=(GEN)sv[3];
  795.     }
  796.     }
  797.   while(k<=n);
  798.   tetpil=avma;return gerepile(av,tetpil,gcopy(u));
  799. }
  800.  
  801. GEN lllgramintindep(x)
  802.      GEN x;
  803. {
  804.   long av=avma,tetpil,lx=lg(x),tx=typ(x),i,j,k,l,p1,n,lim;
  805.   GEN u,B,lam,q,r,h,la,bb,sv;
  806.  
  807.   if(tx!=19) err(lllger1);
  808.   n=lx-1;if(n<=1) return idmat(n);
  809.   if(lg(x[1])!=lx) err(lllger2);
  810.   av=avma;lim=(avma+bot)>>1;
  811.   B=cgetg(lx+1,18);
  812.   B[1]=un;lam=cgetg(lx,19);
  813.   for(j=1;j<=n;j++) lam[j]=lgetg(lx,18);
  814.   for(i=1;i<=n;i++) 
  815.     for(j=1;j<=i;j++)
  816.       {
  817.     u=(GEN)coeff(x,i,j);
  818.     if(typ(u)!=1) err(lllger4);
  819.     for(k=1;k<j;k++)
  820.       u=divii(subii(mulii(B[k+1],u),mulii(coeff(lam,i,k),coeff(lam,j,k))),B[k]);
  821.     if(j<i) coeff(lam,i,j)=(long)u;
  822.     else 
  823.       {
  824.         if(signe(u)) B[i+1]=(long)u;
  825.         else err(lllger5);
  826.       }
  827.     coeff(lam,j,i)=zero;
  828.       }
  829.   k=2;h=idmat(n);
  830.   for(;;)
  831.     {
  832.       if(cmpii(absi(u=shifti(coeff(lam,k,k-1),1)),B[k])>0)
  833.     {
  834.       q=dvmdii(addii(u,B[k]),shifti(B[k],1),&r);
  835.       if(signe(r)<0) q=addsi(-1,q);
  836.       h[k]=lsub(h[k],gmul(q,h[k-1]));
  837.       coeff(lam,k,k-1)=lsubii(coeff(lam,k,k-1),mulii(q,B[k]));
  838.       for(i=1;i<k-1;i++) coeff(lam,k,i)=lsubii(coeff(lam,k,i),mulii(q,coeff(lam,k-1,i)));
  839.     }
  840.       if(cmpii(mulsi(99,mulii(B[k],B[k])),mulsi(100,addii(mulii(B[k-1],B[k+1]),mulii(coeff(lam,k,k-1),coeff(lam,k,k-1)))))>0)
  841.     {
  842.       la=(GEN)coeff(lam,k,k-1);p1=h[k-1];h[k-1]=h[k];h[k]=p1;
  843.       for(j=1;j<=k-2;j++) 
  844.         {
  845.           p1=coeff(lam,k-1,j);coeff(lam,k-1,j)=coeff(lam,k,j);
  846.           coeff(lam,k,j)=p1;
  847.         }
  848.       for(i=k+1;i<=n;i++)
  849.         {
  850.           bb=(GEN)coeff(lam,i,k);
  851.           coeff(lam,i,k)=ldivii(subii(mulii(B[k+1],coeff(lam,i,k-1)),mulii(la,bb)),B[k]);
  852.           coeff(lam,i,k-1)=ldivii(addii(mulii(la,coeff(lam,i,k-1)),mulii(B[k-1],bb)),B[k]);
  853.         }
  854.       B[k]=ldivii(addii(mulii(B[k-1],B[k+1]),mulii(la,la)),B[k]);
  855.       if(k>2) k--;
  856.     }
  857.       else
  858.     {
  859.       for(l=k-2;l>=1;l--)
  860.         {
  861.           if(cmpii(absi(u=shifti(coeff(lam,k,l),1)),B[l+1])>0)
  862.         {
  863.           q=dvmdii(addii(u,B[l+1]),shifti(B[l+1],1),&r);
  864.           if(signe(r)<0) q=addsi(-1,q);
  865.           h[k]=lsub(h[k],gmul(q,h[l]));
  866.           coeff(lam,k,l)=lsubii(coeff(lam,k,l),mulii(q,B[l+1]));
  867.           for(i=1;i<l;i++) coeff(lam,k,i)=lsubii(coeff(lam,k,i),mulii(q,coeff(lam,l,i)));
  868.         }
  869.         }
  870.       k++;
  871.       if(k>n) {tetpil=avma;return gerepile(av,tetpil,gcopy(h));}
  872.     }
  873.       if(avma<lim)
  874.     {
  875.       tetpil=avma;
  876.       sv=cgetg(4,17);
  877.       sv[1]=lcopy(B);sv[2]=lcopy(h);sv[3]=lcopy(lam);
  878.       sv=gerepile(av,tetpil,sv);
  879.       B=(GEN)sv[1];h=(GEN)sv[2];lam=(GEN)sv[3];
  880.     }
  881.  
  882.     }
  883. }
  884.  
  885. GEN lllgramall(x,all)
  886.      GEN x;
  887.      long all;
  888. {
  889.   long av=avma,tetpil,lx=lg(x),tx=typ(x),i,j,k,l,p1,n,lim,dec;
  890.   GEN u,B,lam,q,r,h,la,bb,p2,p3,p4,y,fl;
  891.  
  892.   if(tx!=19) err(lllger1);
  893.   n=lx-1;if(n<=1) return idmat(n);
  894.   if(lg(x[1])!=lx) err(lllger2);
  895.   av=avma;lim=(avma+bot)>>1;
  896.   B=cgetg(lx+1,18);
  897.   fl=cgetg(lx,17);
  898.   B[1]=un;lam=cgetg(lx,19);
  899.   for(j=1;j<lx;j++) lam[j]=lgetg(lx,18);
  900.   for(i=1;i<lx;i++) 
  901.     for(j=1;j<=i;j++)
  902.       {
  903.     if((j<i)&&(!signe(fl[j]))) coeff(lam,i,j)=coeff(lam,j,i)=zero;
  904.     else
  905.       {
  906.         u=(GEN)coeff(x,i,j);
  907.         if(typ(u)!=1) err(lllger4);
  908.         for(k=1;k<j;k++)
  909.           if(signe(fl[k])) u=divii(subii(mulii(B[k+1],u),mulii(coeff(lam,i,k),coeff(lam,j,k))),B[k]);
  910.         if(j<i) {coeff(lam,i,j)=(long)u;coeff(lam,j,i)=zero;}
  911.         else 
  912.           {
  913.         if(signe(u)) {B[i+1]=(long)u;coeff(lam,i,i)=fl[i]=un;}
  914.         else {B[i+1]=B[i];coeff(lam,i,i)=fl[i]=zero;}
  915.           }
  916.       }
  917.       }
  918.   k=2;h=idmat(n);
  919.   for(;;)
  920.     {
  921.       if(cmpii(absi(u=shifti(coeff(lam,k,k-1),1)),B[k])>0)
  922.     {
  923.       q=dvmdii(addii(u,B[k]),shifti(B[k],1),&r);
  924.       if(signe(r)<0) q=addsi(-1,q);
  925.       h[k]=lsub(h[k],gmul(q,h[k-1]));
  926.       coeff(lam,k,k-1)=lsubii(coeff(lam,k,k-1),mulii(q,B[k]));
  927.       for(i=1;i<k-1;i++) coeff(lam,k,i)=lsubii(coeff(lam,k,i),mulii(q,coeff(lam,k-1,i)));
  928.     }
  929.       if(signe(fl[k-1])&&((cmpii(mulsi(99,mulii(B[k],B[k])),mulsi(100,addii(p3=mulii(B[k-1],B[k+1]),p4=mulii(la=(GEN)coeff(lam,k,k-1),coeff(lam,k,k-1)))))>0)||(!signe(fl[k]))))
  930.     {
  931.       p1=h[k-1];h[k-1]=h[k];h[k]=p1;
  932.       for(j=1;j<=k-2;j++) 
  933.         {
  934.           p1=coeff(lam,k-1,j);coeff(lam,k-1,j)=coeff(lam,k,j);
  935.           coeff(lam,k,j)=p1;
  936.         }
  937.       if(signe(fl[k]))
  938.         {
  939.           for(i=k+1;i<=n;i++)
  940.         {
  941.           bb=(GEN)coeff(lam,i,k);
  942.           coeff(lam,i,k)=ldivii(subii(mulii(B[k+1],coeff(lam,i,k-1)),mulii(la,bb)),B[k]);
  943.           coeff(lam,i,k-1)=ldivii(addii(mulii(la,coeff(lam,i,k-1)),mulii(B[k-1],bb)),B[k]);
  944.         }
  945.           B[k]=ldivii(addii(p3,p4),B[k]);
  946.         }
  947.       else
  948.         {
  949.           if(signe(la))
  950.         {
  951.           p2=(GEN)B[k];p1=ldivii(p4,p2);
  952.           for(i=k+1;i<lx;i++)
  953.             coeff(lam,i,k-1)=ldivii(mulii(la,coeff(lam,i,k-1)),p2);
  954.           for(j=k+1;j<lx-1;j++)
  955.             for(i=j+1;i<lx;i++)
  956.               coeff(lam,i,j)=ldivii(mulii(p1,coeff(lam,i,j)),p2);
  957.           B[k+1]=B[k]=p1;
  958.           for(i=k+2;i<=lx;i++)
  959.             B[i]=ldivii(mulii(p1,B[i]),p2);
  960.         }
  961.           else
  962.         {
  963.           coeff(lam,k,k-1)=zero;
  964.           for(i=k+1;i<lx;i++)
  965.             {coeff(lam,i,k)=coeff(lam,i,k-1);coeff(lam,i,k-1)=zero;}
  966.           B[k]=B[k-1];fl[k]=un;fl[k-1]=zero;
  967.         }
  968.         }
  969.       if(k>2) k--;
  970.     }
  971.       else
  972.     {
  973.       for(l=k-2;l>=1;l--)
  974.         {
  975.           if(cmpii(absi(u=shifti(coeff(lam,k,l),1)),B[l+1])>0)
  976.         {
  977.           q=dvmdii(addii(u,B[l+1]),shifti(B[l+1],1),&r);
  978.           if(signe(r)<0) q=addsi(-1,q);
  979.           h[k]=lsub(h[k],gmul(q,h[l]));
  980.           coeff(lam,k,l)=lsubii(coeff(lam,k,l),mulii(q,B[l+1]));
  981.           for(i=1;i<l;i++) coeff(lam,k,i)=lsubii(coeff(lam,k,i),mulii(q,coeff(lam,l,i)));
  982.         }
  983.         }
  984.       k++;
  985.       if(k>n) 
  986.         {
  987.           for(k=1;(k<=n)&&(!signe(fl[k]));k++);
  988.           tetpil=avma;
  989.           if(!all)
  990.         {
  991.           y=cgetg(3,17);
  992.           p2=cgetg(k,19);for(i=1;i<k;i++) p2[i]=lcopy(h[i]);
  993.           y[1]=(long)p2;p2=cgetg(n-k+2,19);y[2]=(long)p2;
  994.           for(i=k;i<=n;i++) p2[i-k+1]=lcopy(h[i]);
  995.         }
  996.           else
  997.         {
  998.           if(all==1)
  999.             {
  1000.               y=cgetg(k,19);for(i=1;i<k;i++) y[i]=lcopy(h[i]);
  1001.             }
  1002.           else
  1003.             {
  1004.               y=cgetg(n-k+2,19);
  1005.               for(i=k;i<=n;i++) y[i-k+1]=lcopy(h[i]);
  1006.             }
  1007.         }
  1008.           return gerepile(av,tetpil,y);
  1009.         }
  1010.     }
  1011.       if(avma<lim)
  1012.     {
  1013.       tetpil=avma;
  1014.       B=gcopy(B);h=gcopy(h);lam=gcopy(lam);fl=gcopy(fl);
  1015.       dec=lpile(av,tetpil,0)>>2;
  1016.       B+=dec;h+=dec;lam+=dec;fl+=dec;
  1017.     }
  1018.     }
  1019. }
  1020.  
  1021.  
  1022. GEN oldlllgramall0(x,all)
  1023.      GEN x;
  1024.      long all;
  1025. {
  1026.   long av=avma,tetpil,lx=lg(x),tx=typ(x),i,j,k,l,p1,n,lim,dec;
  1027.   GEN u,B,lam,q,r,h,la,p2,p3,p4,y,fl;
  1028.  
  1029.   if(tx!=19) err(lllger1);
  1030.   n=lx-1;if(n<=1) return idmat(n);
  1031.   if(lg(x[1])!=lx) err(lllger2);
  1032.   av=avma;lim=(avma+bot)>>1;
  1033.   B=cgetg(lx+1,18);
  1034.   fl=cgetg(lx,17);
  1035.   B[1]=un;lam=cgetg(lx,19);
  1036.   for(j=1;j<lx;j++) lam[j]=lgetg(lx,18);
  1037.   for(i=1;i<lx;i++) 
  1038.     for(j=1;j<=i;j++)
  1039.       {
  1040.     if((j<i)&&(!signe(fl[j]))) coeff(lam,i,j)=coeff(lam,j,i)=zero;
  1041.     else
  1042.       {
  1043.         u=(GEN)coeff(x,i,j);
  1044.         if(typ(u)!=1) err(lllger4);
  1045.         for(k=1;k<j;k++)
  1046.           if(signe(fl[k])) u=divii(subii(mulii(B[k+1],u),mulii(coeff(lam,i,k),coeff(lam,j,k))),B[k]);
  1047.         if(j<i) {coeff(lam,i,j)=(long)u;coeff(lam,j,i)=zero;}
  1048.         else 
  1049.           {
  1050.         if(signe(u)) {B[i+1]=(long)u;coeff(lam,i,i)=fl[i]=un;}
  1051.         else {B[i+1]=B[i];coeff(lam,i,i)=fl[i]=zero;}
  1052.           }
  1053.       }
  1054.       }
  1055.   output(B);
  1056.   k=2;h=idmat(n);
  1057.   for(;;)
  1058.     {
  1059.       printf(" %ld",k);fflush(stdout);
  1060.       if(signe(fl[k-1])&&(!signe(fl[k])))
  1061.     {
  1062.       if(cmpii(absi(u=shifti(coeff(lam,k,k-1),1)),B[k])>0)
  1063.         {
  1064.           q=dvmdii(addii(u,B[k]),shifti(B[k],1),&r);
  1065.           if(signe(r)<0) q=addsi(-1,q);
  1066.           h[k]=lsub(h[k],gmul(q,h[k-1]));
  1067.           coeff(lam,k,k-1)=lsubii(coeff(lam,k,k-1),mulii(q,B[k]));
  1068.           for(i=1;i<k-1;i++) coeff(lam,k,i)=lsubii(coeff(lam,k,i),mulii(q,coeff(lam,k-1,i)));
  1069.         }
  1070.       la=(GEN)coeff(lam,k,k-1);
  1071.       p3=mulii(B[k-1],B[k+1]);p4=mulii(la,la);
  1072.       p1=h[k-1];h[k-1]=h[k];h[k]=p1;
  1073.       for(j=1;j<=k-2;j++) 
  1074.         {
  1075.           p1=coeff(lam,k-1,j);coeff(lam,k-1,j)=coeff(lam,k,j);
  1076.           coeff(lam,k,j)=p1;
  1077.         }
  1078.       if(signe(la))
  1079.         {
  1080.           p2=(GEN)B[k];p1=ldivii(p4,p2);
  1081.           for(i=k+1;i<lx;i++)
  1082.         coeff(lam,i,k-1)=ldivii(mulii(la,coeff(lam,i,k-1)),p2);
  1083.           for(j=k+1;j<lx-1;j++)
  1084.         for(i=j+1;i<lx;i++)
  1085.           coeff(lam,i,j)=ldivii(mulii(p1,coeff(lam,i,j)),p2);
  1086.           B[k+1]=B[k]=p1;
  1087.           for(i=k+2;i<=lx;i++)
  1088.         B[i]=ldivii(mulii(p1,B[i]),p2);
  1089.         }
  1090.       else
  1091.         {
  1092.           for(i=k+1;i<lx;i++)
  1093.         {coeff(lam,i,k)=coeff(lam,i,k-1);coeff(lam,i,k-1)=zero;}
  1094.           B[k]=B[k-1];fl[k]=un;fl[k-1]=zero;
  1095.         }
  1096.       if(k>2) k--;
  1097.     }
  1098.       else
  1099.     {
  1100.       for(l=k-1;l>=1;l--)
  1101.         {
  1102.           if(cmpii(absi(u=shifti(coeff(lam,k,l),1)),B[l+1])>0)
  1103.         {
  1104.           q=dvmdii(addii(u,B[l+1]),shifti(B[l+1],1),&r);
  1105.           if(signe(r)<0) q=addsi(-1,q);
  1106.           h[k]=lsub(h[k],gmul(q,h[l]));
  1107.           coeff(lam,k,l)=lsubii(coeff(lam,k,l),mulii(q,B[l+1]));
  1108.           for(i=1;i<l;i++) coeff(lam,k,i)=lsubii(coeff(lam,k,i),mulii(q,coeff(lam,l,i)));
  1109.         }
  1110.         }
  1111.       k++;
  1112.       if(k>n) 
  1113.         {
  1114.           for(k=1;(k<=n)&&(!signe(fl[k]));k++);
  1115.           tetpil=avma;
  1116.           if(!all)
  1117.         {
  1118.           y=cgetg(3,17);
  1119.           p2=cgetg(k,19);for(i=1;i<k;i++) p2[i]=lcopy(h[i]);
  1120.           y[1]=(long)p2;p2=cgetg(n-k+2,19);y[2]=(long)p2;
  1121.           for(i=k;i<=n;i++) p2[i-k+1]=lcopy(h[i]);
  1122.         }
  1123.           else
  1124.         {
  1125.           if(all==1)
  1126.             {
  1127.               y=cgetg(k,19);for(i=1;i<k;i++) y[i]=lcopy(h[i]);
  1128.             }
  1129.           else
  1130.             {
  1131.               y=cgetg(n-k+2,19);
  1132.               for(i=k;i<=n;i++) y[i-k+1]=lcopy(h[i]);
  1133.             }
  1134.         }
  1135.           return gerepile(av,tetpil,y);
  1136.         }
  1137.     }
  1138.       if(avma<lim)
  1139.     {
  1140.       tetpil=avma;
  1141.       B=gcopy(B);h=gcopy(h);lam=gcopy(lam);fl=gcopy(fl);
  1142.       dec=lpile(av,tetpil,0)>>2;
  1143.       B+=dec;h+=dec;lam+=dec;fl+=dec;
  1144.     }
  1145.     }
  1146. }
  1147.  
  1148. GEN lllall0(x,all)
  1149.      GEN x;
  1150.      long all;
  1151. {
  1152.   long av=avma,tetpil,lx=lg(x),tx=typ(x),i,j,k,l,p1,n,lim,dec,kmax;
  1153.   GEN u,B,lam,q,r,h,la,p2,p3,p4,y,fl;
  1154.  
  1155.   if(tx!=19) err(lllger1);
  1156.   n=lx-1;
  1157.   if(!n) 
  1158.     {
  1159.       if(all) return cgetg(1,19);
  1160.       else {y=cgetg(3,17);y[1]=lgetg(1,19);y[2]=lgetg(1,19);return y;}
  1161.     }
  1162.   if(n==1) 
  1163.     {
  1164.       if(gcmp0(x[1]))
  1165.     {
  1166.       if(!all) 
  1167.         {y=cgetg(3,17);y[1]=(long)idmat(1);y[2]=lgetg(1,19);return y;}
  1168.       return (all==1)?idmat(1):cgetg(1,19);
  1169.     }
  1170.       else
  1171.     {
  1172.       if(!all) {y=cgetg(3,17);y[1]=lgetg(1,19);y[2]=lcopy(x);return y;}
  1173.       return (all==1)?cgetg(1,19):gcopy(x);
  1174.     }
  1175.     }
  1176.   av=avma;x=gcopy(x);lim=(avma+bot)>>1;
  1177.   B=cgetg(lx+1,18);for(i=1;i<=lx;i++) B[i]=zero;
  1178.   fl=cgetg(lx,17);for(i=1;i<lx;i++) fl[i]=zero;
  1179.   B[1]=un;lam=idmat(n);
  1180.   u=gscal(x[1],x[1]);
  1181.   if(signe(u)) {B[2]=(long)u;coeff(lam,1,1)=fl[1]=un;}
  1182.   else {B[2]=un;coeff(lam,1,1)=fl[1]=zero;}
  1183.   k=2;h=idmat(n);kmax=1;
  1184.   for(;;)
  1185.     {
  1186.       if(k>kmax)
  1187.     {
  1188.       kmax=k;
  1189.       for(j=1;j<=k;j++)
  1190.         {
  1191.           if((j<k)&&(!signe(fl[j]))) coeff(lam,k,j)=zero;
  1192.           else
  1193.         {
  1194.           u=gscal(x[k],x[j]);
  1195.           if(typ(u)!=1) err(lllger4);
  1196.           for(i=1;i<j;i++)
  1197.             if(signe(fl[i])) u=divii(subii(mulii(B[i+1],u),mulii(coeff(lam,k,i),coeff(lam,j,i))),B[i]);
  1198.           if(j<k) coeff(lam,k,j)=(long)u;
  1199.           else 
  1200.             {
  1201.               if(signe(u)) {B[k+1]=(long)u;coeff(lam,k,k)=fl[k]=un;}
  1202.               else {B[k+1]=B[k];coeff(lam,k,k)=fl[k]=zero;}
  1203.             }
  1204.         }
  1205.         }
  1206.     }
  1207.       if(signe(fl[k-1])&&(!signe(fl[k])))
  1208.     {
  1209.       if(cmpii(absi(u=shifti(coeff(lam,k,k-1),1)),B[k])>0)
  1210.         {
  1211.           q=dvmdii(addii(u,B[k]),shifti(B[k],1),&r);
  1212.           if(signe(r)<0) q=addsi(-1,q);
  1213.           h[k]=lsub(h[k],gmul(q,h[k-1]));
  1214.           x[k]=lsub(x[k],gmul(q,x[k-1]));
  1215.           coeff(lam,k,k-1)=lsubii(coeff(lam,k,k-1),mulii(q,B[k]));
  1216.           for(i=1;i<k-1;i++) coeff(lam,k,i)=lsubii(coeff(lam,k,i),mulii(q,coeff(lam,k-1,i)));
  1217.         }
  1218.       la=(GEN)coeff(lam,k,k-1);
  1219.       p3=mulii(B[k-1],B[k+1]);p4=mulii(la,la);
  1220.       p1=h[k-1];h[k-1]=h[k];h[k]=p1;
  1221.       p1=x[k-1];x[k-1]=x[k];x[k]=p1;
  1222.       for(j=1;j<=k-2;j++) 
  1223.         {
  1224.           p1=coeff(lam,k-1,j);coeff(lam,k-1,j)=coeff(lam,k,j);
  1225.           coeff(lam,k,j)=p1;
  1226.         }
  1227.       if(signe(la))
  1228.         {
  1229.           p2=(GEN)B[k];p1=ldivii(p4,p2);
  1230.           for(i=k+1;i<=kmax;i++)
  1231.         coeff(lam,i,k-1)=ldivii(mulii(la,coeff(lam,i,k-1)),p2);
  1232.           for(j=k+1;j<kmax;j++)
  1233.         for(i=j+1;i<=kmax;i++)
  1234.           coeff(lam,i,j)=ldivii(mulii(p1,coeff(lam,i,j)),p2);
  1235.           B[k+1]=B[k]=p1;
  1236.           for(i=k+2;i<=kmax+1;i++)
  1237.         B[i]=ldivii(mulii(p1,B[i]),p2);
  1238.         }
  1239.       else 
  1240.         {
  1241.           for(i=k+1;i<=kmax;i++)
  1242.         {coeff(lam,i,k)=coeff(lam,i,k-1);coeff(lam,i,k-1)=zero;}
  1243.           B[k]=B[k-1];fl[k]=un;fl[k-1]=zero;
  1244.         }
  1245.       if(k>2) k--;
  1246.     }
  1247.       else
  1248.     {
  1249.       for(l=k-1;l>=1;l--)
  1250.         {
  1251.           if(cmpii(absi(u=shifti(coeff(lam,k,l),1)),B[l+1])>0)
  1252.         {
  1253.           q=dvmdii(addii(u,B[l+1]),shifti(B[l+1],1),&r);
  1254.           if(signe(r)<0) q=addsi(-1,q);
  1255.           h[k]=lsub(h[k],gmul(q,h[l]));
  1256.           x[k]=lsub(x[k],gmul(q,x[l]));
  1257.           coeff(lam,k,l)=lsubii(coeff(lam,k,l),mulii(q,B[l+1]));
  1258.           for(i=1;i<l;i++) coeff(lam,k,i)=lsubii(coeff(lam,k,i),mulii(q,coeff(lam,l,i)));
  1259.         }
  1260.         }
  1261.       k++;
  1262.       if(k>n) 
  1263.         {
  1264.           for(k=1;(k<=n)&&(!signe(fl[k]));k++);
  1265.           tetpil=avma;
  1266.           if(!all)
  1267.         {
  1268.           y=cgetg(3,17);
  1269.           p2=cgetg(k,19);for(i=1;i<k;i++) p2[i]=lcopy(h[i]);
  1270.           y[1]=(long)p2;p2=cgetg(n-k+2,19);y[2]=(long)p2;
  1271.           for(i=k;i<=n;i++) p2[i-k+1]=lcopy(h[i]);
  1272.         }
  1273.           else
  1274.         {
  1275.           if(all==1)
  1276.             {
  1277.               y=cgetg(k,19);for(i=1;i<k;i++) y[i]=lcopy(h[i]);
  1278.             }
  1279.           else
  1280.             {
  1281.               y=cgetg(n-k+2,19);
  1282.               for(i=k;i<=n;i++) y[i-k+1]=lcopy(h[i]);
  1283.             }
  1284.         }
  1285.           return gerepile(av,tetpil,y);
  1286.         }
  1287.     }
  1288.       if(avma<lim)
  1289.     {
  1290.       tetpil=avma;
  1291.       B=gcopy(B);h=gcopy(h);lam=gcopy(lam);fl=gcopy(fl);x=gcopy(x);
  1292.       dec=lpile(av,tetpil,0)>>2;
  1293.       B+=dec;h+=dec;lam+=dec;fl+=dec;x+=dec;
  1294.     }
  1295.     }
  1296. }
  1297.  
  1298. GEN newlllall0(x,all)
  1299.      GEN x;
  1300.      long all;
  1301. {
  1302.   long av=avma,tetpil,lx=lg(x),tx=typ(x),i,j,k,l,p1,n,lim,dec,kmax;
  1303.   GEN u,B,lam,q,r,h,la,p2,y,fl,a,b,c,d;
  1304.  
  1305.   if(tx!=19) err(lllger1);
  1306.   n=lx-1;
  1307.   if(!n) 
  1308.     {
  1309.       if(all) return cgetg(1,19);
  1310.       else {y=cgetg(3,17);y[1]=lgetg(1,19);y[2]=lgetg(1,19);return y;}
  1311.     }
  1312.   if(n==1) 
  1313.     {
  1314.       if(gcmp0(x[1]))
  1315.     {
  1316.       if(!all) 
  1317.         {y=cgetg(3,17);y[1]=(long)idmat(1);y[2]=lgetg(1,19);return y;}
  1318.       return (all==1)?idmat(1):cgetg(1,19);
  1319.     }
  1320.       else
  1321.     {
  1322.       if(!all) {y=cgetg(3,17);y[1]=lgetg(1,19);y[2]=lcopy(x);return y;}
  1323.       return (all==1)?cgetg(1,19):gcopy(x);
  1324.     }
  1325.     }
  1326.   av=avma;x=gcopy(x);lim=(avma+bot)>>1;
  1327.   B=cgetg(lx+1,18);
  1328.   fl=cgetg(lx,17);
  1329.   B[1]=un;lam=cgetg(lx,19);
  1330.   for(j=1;j<lx;j++) lam[j]=lgetg(lx,18);
  1331.   u=gscal(x[1],x[1]);
  1332.   if(signe(u)) {B[2]=(long)u;coeff(lam,1,1)=fl[1]=un;}
  1333.   else {B[2]=un;coeff(lam,1,1)=fl[1]=zero;}
  1334.   k=2;h=idmat(n);kmax=1;
  1335.   for(;;)
  1336.     {
  1337.       if(k>kmax)
  1338.     {
  1339.       kmax=k;
  1340.       for(j=1;j<=k;j++)
  1341.         {
  1342.           if((j<k)&&(!signe(fl[j]))) coeff(lam,k,j)=coeff(lam,j,k)=zero;
  1343.           else
  1344.         {
  1345.           u=gscal(x[k],x[j]);
  1346.           if(typ(u)!=1) err(lllger4);
  1347.           for(i=1;i<j;i++)
  1348.             if(signe(fl[i])) u=divii(subii(mulii(B[i+1],u),mulii(coeff(lam,k,i),coeff(lam,j,i))),B[i]);
  1349.           if(j<k) {coeff(lam,k,j)=(long)u;coeff(lam,j,k)=zero;}
  1350.           else 
  1351.             {
  1352.               if(signe(u)) {B[k+1]=(long)u;coeff(lam,k,k)=fl[k]=un;}
  1353.               else {B[k+1]=B[k];coeff(lam,k,k)=fl[k]=zero;}
  1354.             }
  1355.         }
  1356.         }
  1357.     }
  1358.       printf(" %ld",k);fflush(stdout);
  1359.       if(signe(fl[k-1])&&(!signe(fl[k])))
  1360.     {
  1361. /*      printf("\nswap"); */
  1362.       la=(GEN)coeff(lam,k,k-1);
  1363.       p2=bezout(la,B[k],&a,&b);d=negi(divii(la,p2));c=divii(B[k],p2);
  1364.       p1=ladd(gmul(a,h[k]),gmul(b,h[k-1]));
  1365.       h[k-1]=ladd(gmul(c,h[k]),gmul(d,h[k-1]));h[k]=p1;
  1366.       p1=ladd(gmul(a,x[k]),gmul(b,x[k-1]));
  1367.       x[k-1]=ladd(gmul(c,x[k]),gmul(d,x[k-1]));x[k]=p1;
  1368.       for(j=1;j<=k-2;j++) 
  1369.         {
  1370.           p1=laddii(mulii(a,coeff(lam,k,j)),mulii(b,coeff(lam,k-1,j)));
  1371.           coeff(lam,k-1,j)=laddii(mulii(c,coeff(lam,k,j)),mulii(d,coeff(lam,k-1,j)));
  1372.           coeff(lam,k,j)=p1;
  1373.         }
  1374.       p2=mulii(c,c);
  1375.       for(i=k+1;i<=kmax;i++)
  1376.         {coeff(lam,i,k)=ldivii(coeff(lam,i,k-1),p2);coeff(lam,i,k-1)=zero;}
  1377.       for(j=k+1;j<kmax;j++)
  1378.         for(i=j+1;i<=kmax;i++)
  1379.           coeff(lam,i,j)=ldivii(coeff(lam,i,j),p2);
  1380.       B[k+1]=ldivii(B[k+1],p2);
  1381. /*      printf("\nd_k=");output(gdiv(B[k],p2));  */
  1382.       B[k]=B[k-1];fl[k]=un;fl[k-1]=zero;
  1383.       for(i=k+2;i<=kmax+1;i++)
  1384.         B[i]=ldivii(B[i],p2);
  1385.       if(k>2) k--;
  1386.     }
  1387.       else
  1388.     {
  1389.       for(l=k-1;l>=1;l--)
  1390.         {
  1391.           if(cmpii(absi(u=shifti(coeff(lam,k,l),1)),B[l+1])>0)
  1392.         {
  1393. /*          printf("\nk k-L"); */
  1394.           q=dvmdii(addii(u,B[l+1]),shifti(B[l+1],1),&r);
  1395.           if(signe(r)<0) q=addsi(-1,q);
  1396.           h[k]=lsub(h[k],gmul(q,h[l]));
  1397.           x[k]=lsub(x[k],gmul(q,x[l]));
  1398.           coeff(lam,k,l)=lsubii(coeff(lam,k,l),mulii(q,B[l+1]));
  1399.           for(i=1;i<l;i++) coeff(lam,k,i)=lsubii(coeff(lam,k,i),mulii(q,coeff(lam,l,i)));
  1400.         }
  1401.         }
  1402.       k++;
  1403.       if(k>n) 
  1404.         {
  1405.           for(k=1;(k<=n)&&(!signe(fl[k]));k++);
  1406.           tetpil=avma;
  1407.           if(!all)
  1408.         {
  1409.           y=cgetg(3,17);
  1410.           p2=cgetg(k,19);for(i=1;i<k;i++) p2[i]=lcopy(h[i]);
  1411.           y[1]=(long)p2;p2=cgetg(n-k+2,19);y[2]=(long)p2;
  1412.           for(i=k;i<=n;i++) p2[i-k+1]=lcopy(h[i]);
  1413.         }
  1414.           else
  1415.         {
  1416.           if(all==1)
  1417.             {
  1418.               y=cgetg(k,19);for(i=1;i<k;i++) y[i]=lcopy(h[i]);
  1419.             }
  1420.           else
  1421.             {
  1422.               y=cgetg(n-k+2,19);
  1423.               for(i=k;i<=n;i++) y[i-k+1]=lcopy(h[i]);
  1424.             }
  1425.         }
  1426.           return gerepile(av,tetpil,y);
  1427.         }
  1428.     }
  1429.       if(avma<lim)
  1430.     {
  1431.       tetpil=avma;
  1432.       B=gcopy(B);h=gcopy(h);lam=gcopy(lam);fl=gcopy(fl);x=gcopy(x);
  1433.       dec=lpile(av,tetpil,0)>>2;
  1434.       B+=dec;h+=dec;lam+=dec;fl+=dec;x+=dec;
  1435.     }
  1436.     }
  1437. }
  1438.  
  1439.  
  1440. GEN lllgramint(x)
  1441.      GEN x;
  1442. {
  1443.   return lllgramall(x,2);
  1444. }
  1445.  
  1446. GEN lllgramkerim(x)
  1447.      GEN x;
  1448. {
  1449.   return lllgramall(x,0);
  1450. }
  1451.  
  1452. /********************************************************************/
  1453. /********************************************************************/
  1454. /**                                                                **/
  1455. /**               DEPENDANCE LINEAIRE ET ALGEBRIQUE                **/
  1456. /**                                                                **/
  1457. /********************************************************************/
  1458. /********************************************************************/
  1459.  
  1460. GEN    lindep4(x,bit,prec)    /* C.B. 27/11/91 */
  1461.  
  1462.      GEN  x;
  1463.      long bit,prec;
  1464.      
  1465. {
  1466.   long lx=lg(x),cpt=0;
  1467.   long av0,av,tetpil,lim,l,i,j,k,n,pro,mu1,mu2,e;
  1468.   GEN  G0,G1,G,D,mu,u,B,BB,BK,p,q,r,cst,unreel,sv;
  1469.  
  1470.  
  1471.  
  1472.   av0=avma;
  1473.   cst=gdivgs(stoi(95),100); /* LLL proposent 0.75 */
  1474.   lim=(avma+bot)>>1;
  1475.   if (prec)
  1476.     {
  1477.       affsr(1,unreel=cgetr(prec+1));
  1478.       cst=gmul(cst,unreel);
  1479.     }
  1480.   G0=idmat(n=lx-1);
  1481.   for(i=1;i<=n;i++)
  1482.     {
  1483.       for(j=i+1;j<=n;j++)
  1484.     coeff(G0,i,j)=coeff(G0,j,i)=lshift(greal(gmul(x[i],x[j])),bit);
  1485.       coeff(G0,i,i)=lshift(gnorm(x[i]),bit);
  1486.       if(i>1)
  1487.     coeff(G0,i,i)=laddgs(coeff(G0,i,i),1);
  1488.     }
  1489.   sor(G0,'e',12,0);
  1490.   av=avma;
  1491.  
  1492.  
  1493.   mu=idmat(n=lx-1);
  1494.   B=cgetg(lx,17);
  1495.   B[1]=lnorm(x[1]);
  1496.   B[2]=limag(gmul(x[1],gconj(x[2])));
  1497.   for(j=2;j<lx;j++)
  1498.     coeff(mu,j,1)=ldiv(greal(gmul(x[1],gconj(x[j]))),B[1]);
  1499.   if(!gcmp0(B[2]))
  1500.      {
  1501.        for(j=3;j<lx;j++)
  1502.      coeff(mu,j,2)=ldiv(gimag(gmul(x[1],gconj(x[j]))),B[2]);
  1503.        B[2]=lshift(gdiv(gmul(B[2],B[2]),B[1]),bit);
  1504.      }
  1505.  
  1506.   else B[2 ]=(long)unreel;
  1507.   B[1]=lshift(B[1],bit);
  1508.   for(j=3;j<=n;j++) B[j]=(long)unreel;
  1509.  
  1510.   D=idmat(n);
  1511.   for(i=1;i<=n;i++) coeff(D,i,i)=B[i];
  1512.   G1=gmul(gmul(mu,D),gtrans(mu));
  1513.   sor(G1,'e',12,0);
  1514.  
  1515.  
  1516.   sor(mu,'e',12,0);
  1517.   sor(sqred3(G1),'e',12,0);
  1518. /*-- LLL : ---------------------------------------*/
  1519.  
  1520.   u=idmat(n);
  1521.   k=2;cpt=0;
  1522.   do
  1523.     {
  1524. /*      printf("k=%d  cpt=%d\n",k,cpt);fflush(stdout); */
  1525.       if(!gcmp0(r=grndtoi(coeff(mu,k,k-1),&e)))
  1526.         {
  1527.           u[k]=lsub(u[k],gmul(r,u[k-1]));
  1528.           for(j=1;j<k-1;j++)
  1529.             coeff(mu,k,j)=lsub(coeff(mu,k,j),gmul(r,coeff(mu,k-1,j)));
  1530.       mu1=coeff(mu,k,k-1)=lsub(coeff(mu,k,k-1),r);
  1531.       affrr(mu1,pro=lgetr(prec));
  1532.       mu1=coeff(mu,k,k-1)=pro;
  1533.         }
  1534.       else mu1=coeff(mu,k,k-1);
  1535.       q=gmul(B[k-1],gsub(cst,mu2=lsqr(mu1)));
  1536.       if(gcmp(q,B[k])>0)
  1537.         {
  1538. /********** ici on permute bk et b(k-1) ***************************/
  1539.  
  1540.           BB=gadd(B[k],gmul(B[k-1],mu2));
  1541.           coeff(mu,k,k-1)=ldiv(gmul(mu1,B[k-1]),BB);
  1542.           B[k]=lmul(B[k-1],BK=gdiv(B[k],BB));
  1543.           B[k-1]=(long)BB;
  1544.           pro=u[k];u[k]=u[k-1];u[k-1]=pro;
  1545.           for(j=1;j<=k-2;j++)
  1546.             {
  1547.               pro=coeff(mu,k,j);coeff(mu,k,j)=coeff(mu,k-1,j);
  1548.               coeff(mu,k-1,j)=pro;
  1549.             }
  1550.           for(i=k+1;i<=n;i++)
  1551.             {
  1552.               p=(GEN)coeff(mu,i,k);
  1553.           coeff(mu,i,k)=lsub(coeff(mu,i,k-1),gmul(mu1,p));
  1554.               coeff(mu,i,k-1)=ladd(gmul(BK,p),gmul(coeff(mu,k,k-1),coeff(mu,i,k-1)));
  1555.             }
  1556.           if(k>2) {k--;cpt++;}
  1557.         }
  1558.       else
  1559.         {
  1560. /********** ici bk et b(k-1) OK on rend propre *******************/
  1561.           for(l=k-2;l;l--)
  1562.         if(!gcmp0(r=grndtoi(coeff(mu,k,l),&e)))
  1563.           {
  1564.         u[k]=lsub(u[k],gmul(r,u[l]));
  1565.         for(j=1;j<l;j++)
  1566.           coeff(mu,k,j)=lsub(coeff(mu,k,j),gmul(r,coeff(mu,l,j)));
  1567.         
  1568.         coeff(mu,k,l)=lsub(coeff(mu,k,l),r);
  1569.           }
  1570.       G=gmul(gmul(gtrans(u),G0),u);
  1571.   printf("sqred1 : cpt = %d\n",cpt);
  1572.       mu=gtrans(sqred1(G));
  1573.           k++;cpt++;
  1574. /******************************************************************/
  1575.         }
  1576.       if(avma<lim)
  1577.     {
  1578.       tetpil=avma;
  1579.       sv=cgetg(4,17);
  1580.       sv[1]=lcopy(B);sv[2]=lcopy(u);sv[3]=lcopy(mu);
  1581.       sv=gerepile(av,tetpil,sv);
  1582.       B=(GEN)sv[1];u=(GEN)sv[2];mu=(GEN)sv[3];
  1583.     }
  1584.     }
  1585.   while(k<=n);
  1586.   printf("\n *** %d etapes dans LLL\n",cpt);
  1587.   tetpil=avma;return gerepile(av,tetpil,gcopy(u[1]));
  1588. }
  1589.  
  1590. GEN    lindep3(x,bit,prec)    /* C.B. 13/11/91 */
  1591.  
  1592.      GEN  x;
  1593.      long bit,prec;
  1594.      
  1595. {
  1596.   long lx=lg(x),cpt=0;
  1597.   long av0,av,tetpil,lim,l,i,j,k,n,pro,mu1,mu2,e;
  1598.   GEN  mu,u,B,BB,BK,p,q,r,cst,unreel,sv;
  1599.  
  1600.  
  1601.  
  1602.   av0=avma;
  1603.   cst=gdivgs(stoi(95),100); /* LLL proposent 0.75 */
  1604.   lim=(avma+bot)>>1;
  1605.   if (prec)
  1606.     {
  1607.       affsr(1,unreel=cgetr(prec+1));
  1608.       cst=gmul(cst,unreel);
  1609.     }
  1610.   av=avma;
  1611.  
  1612.   mu=idmat(n=lx-1);
  1613.   B=cgetg(lx,17);
  1614.   B[1]=lnorm(x[1]);
  1615.   B[2]=limag(gmul(x[1],gconj(x[2])));
  1616.   for(j=2;j<lx;j++)
  1617.     coeff(mu,j,1)=ldiv(greal(gmul(x[1],gconj(x[j]))),B[1]);
  1618.   if(!gcmp0(B[2]))
  1619.      {
  1620.        for(j=3;j<lx;j++)
  1621.      coeff(mu,j,2)=ldiv(gimag(gmul(x[1],gconj(x[j]))),B[2]);
  1622.        B[2]=lshift(gdiv(gmul(B[2],B[2]),B[1]),bit);
  1623.      }
  1624.  
  1625.   else B[2 ]=(long)unreel;
  1626.   B[1]=lshift(B[1],bit);
  1627.   for(j=3;j<=n;j++) B[j]=(long)unreel;
  1628. /*-- LLL : ---------------------------------------*/
  1629.  
  1630.   u=idmat(n);
  1631.   k=2;cpt=0;
  1632.   do
  1633.     {
  1634.  
  1635. /*      printf("k=%d  cpt=%d\n",k,cpt);fflush(stdout); */
  1636.       if(!gcmp0(r=grndtoi(coeff(mu,k,k-1),&e)))
  1637.         {
  1638.           u[k]=lsub(u[k],gmul(r,u[k-1]));
  1639.           for(j=1;j<k-1;j++)
  1640.             coeff(mu,k,j)=lsub(coeff(mu,k,j),gmul(r,coeff(mu,k-1,j)));
  1641.       mu1=coeff(mu,k,k-1)=lsub(coeff(mu,k,k-1),r);
  1642.       affrr(mu1,pro=lgetr(prec));
  1643.       mu1=coeff(mu,k,k-1)=pro;
  1644.         }
  1645.       else mu1=coeff(mu,k,k-1);
  1646.       q=gmul(B[k-1],gsub(cst,mu2=lsqr(mu1)));
  1647.       if(gcmp(q,B[k])>0)
  1648.         {
  1649.           BB=gadd(B[k],gmul(B[k-1],mu2));
  1650.           coeff(mu,k,k-1)=ldiv(gmul(mu1,B[k-1]),BB);
  1651.           B[k]=lmul(B[k-1],BK=gdiv(B[k],BB));
  1652.           B[k-1]=(long)BB;
  1653.           pro=u[k];u[k]=u[k-1];u[k-1]=pro;
  1654.           for(j=1;j<=k-2;j++)
  1655.             {
  1656.               pro=coeff(mu,k,j);coeff(mu,k,j)=coeff(mu,k-1,j);
  1657.               coeff(mu,k-1,j)=pro;
  1658.             }
  1659.           for(i=k+1;i<=n;i++)
  1660.             {
  1661.               p=(GEN)coeff(mu,i,k);
  1662.           coeff(mu,i,k)=lsub(coeff(mu,i,k-1),gmul(mu1,p));
  1663.               coeff(mu,i,k-1)=ladd(gmul(BK,p),gmul(coeff(mu,k,k-1),coeff(mu,i,k-1)));
  1664.             }
  1665.           if(k>2) {k--;cpt++;}
  1666.         }
  1667.       else
  1668.         {
  1669.           for(l=k-2;l;l--)
  1670.             {
  1671.               if(!gcmp0(r=grndtoi(coeff(mu,k,l),&e)))
  1672.                 {
  1673. /*          printf("k, l, e, r: %d, %d, %d, ",k,l,e);output(r); */
  1674.                   u[k]=lsub(u[k],gmul(r,u[l]));
  1675.                   for(j=1;j<l;j++)
  1676.                     coeff(mu,k,j)=lsub(coeff(mu,k,j),gmul(r,coeff(mu,l,j)));
  1677.  
  1678.           pro=coeff(mu,k,l)=lsub(coeff(mu,k,l),r);
  1679.           affrr(pro,coeff(mu,k,l)=lgetr(prec));
  1680.                 }
  1681.             }
  1682.           k++;cpt++;
  1683.         }
  1684.       if(avma<lim)
  1685.     {
  1686.       tetpil=avma;
  1687.       sv=cgetg(4,17);
  1688.       sv[1]=lcopy(B);sv[2]=lcopy(u);sv[3]=lcopy(mu);
  1689.       sv=gerepile(av,tetpil,sv);
  1690.       B=(GEN)sv[1];u=(GEN)sv[2];mu=(GEN)sv[3];
  1691.     }
  1692.     }
  1693.   while(k<=n);
  1694.   printf("\n *** %d etapes dans LLL\n",cpt);
  1695.   tetpil=avma;return gerepile(av,tetpil,gcopy(u[1]));
  1696. }
  1697.  
  1698.  
  1699. GEN    lindep2(x,bit,prec)
  1700.      
  1701.      GEN  x;
  1702.      long bit,prec;
  1703.      
  1704. {
  1705.   long  tx=typ(x),lx=lg(x),ly,i,j,flag,av,tetpil,e;
  1706.   GEN   y,p1,p2,p3,p4,p5;
  1707.   
  1708.   if((tx<17)||(tx>18)) err(linder1);
  1709.   av=avma;p1=greal(x);p2=gimag(x);
  1710.   ly=(flag=!gcmp0(p2)) ? lx+2 : lx+1;
  1711.   p4=cgetg(lx,19);
  1712.   for(i=1;i<lx;i++)
  1713.   {
  1714.     p5=cgetg(ly,18);p4[i]=(long)p5;
  1715.     for(j=1;j<lx;j++) p5[j]=(i==j) ? un : zero;
  1716.     p5[lx]=lcvtoi(gshift(p1[i],bit),&e);
  1717.     if(flag) p5[lx+1]=lcvtoi(gshift(p2[i],bit),&e);
  1718.   }
  1719.   p5=gmul(p4,lllint(p4));p3=(GEN)p5[1];
  1720.   tetpil=avma;y=cgetg(lx,17);
  1721.   for(i=1;i<lx;i++) y[i]=lcopy(p3[i]);
  1722.   y=gerepile(av,tetpil,y);
  1723.   return y;
  1724. }
  1725.  
  1726. #define quazero(x) (gcmp0(x)||((typ(x)==2)&&(expo(x)< -32*(prec-2)+2*n)))
  1727.  
  1728. GEN    lindep(x,prec)
  1729.      
  1730.      GEN  x;
  1731.      long prec;
  1732.      
  1733. {
  1734.   GEN b[50],be[50],bn[50],m[50][50],y,c1,c2,c3,px,py,pxy;
  1735.   GEN p1,p2,p3,p4,p5,p6,p7,r,f,em;
  1736.   long av,av1,tetpil,lx=lg(x),tx=typ(x),n,i,j,fl,i1;
  1737.   long qzer[50];
  1738.   
  1739.   if((tx<17)||(tx>18)) err(linder1);
  1740.   if(lx>=50) err(linder2);
  1741.   av=avma;n=lx-1;p1=greal(x);p2=gimag(x);
  1742.   px=gscal(p1,p1);py=gscal(p2,p2);
  1743.   pxy=gscal(p1,p2);p3=mpsub(mpmul(px,py),mpmul(pxy,pxy));
  1744.   for(i=1;i<=n;i++)
  1745.   {
  1746.     be[i]=cgetg(lx,18);
  1747.     for(j=1;j<=n;j++) be[i][j]=lgetr(prec+1);
  1748.   }
  1749.   for(j=1;j<=n;j++) bn[j]=cgetr(prec+1);
  1750.   for(i=1;i<=n;i++)
  1751.   {
  1752.     for(j=1;j<i;j++)
  1753.       m[i][j]=cgetr(prec+1);
  1754.   }
  1755.   if(quazero(p1)) {p1=p2;px=py;fl=1;} else fl=quazero(p3);
  1756.   for(i=1;i<=n;i++)
  1757.   {
  1758.     b[i]=cgetg(lx,18);
  1759.     for(j=1;j<=n;j++)
  1760.       b[i][j]=(i==j) ? un : zero;
  1761.   }
  1762.   av1=avma;
  1763.   for(i=1;i<=n;i++)
  1764.   {
  1765.     if(fl) p4=gmul(gdiv(gscal(b[i],p1),px),p1);
  1766.     else
  1767.     {
  1768.       p4=gscal(b[i],p1);p5=gscal(b[i],p2);
  1769.       p6=gdiv(mpsub(mpmul(py,p4),mpmul(pxy,p5)),p3);
  1770.       p7=gdiv(mpsub(mpmul(px,p5),mpmul(pxy,p4)),p3);
  1771.       p4=gadd(gmul(p6,p1),gmul(p7,p2));
  1772.     }
  1773.     if(tx==18) p4=gsub(b[i],p4);
  1774.     else p4=gsub(b[i],gtrans(p4));
  1775.     for(j=1;j<i;j++)
  1776.       if(!qzer[j])
  1777.       {
  1778.         gdivz(gscal(b[i],be[j]),bn[j],m[i][j]);
  1779.         p4=gsub(p4,gmul(m[i][j],be[j]));
  1780.       }
  1781.       else affrr(bn[j],m[i][j]);
  1782.     gaffect(p4,be[i]);
  1783.     affrr(gscal(be[i],be[i]),bn[i]);
  1784.     qzer[i]=quazero(bn[i]);avma=av1;
  1785.   }
  1786.   while(qzer[n])
  1787.   {
  1788.     em=bn[1];j=1;av1=avma;
  1789.     for(i=2;i<n;i++)
  1790.     {
  1791.       p3=shiftr(bn[i],i);
  1792.       if(cmprr(p3,em)>0)
  1793.       {
  1794.         em=p3;j=i;
  1795.       }
  1796.     }
  1797.     avma=av1;i=j;i1=i+1;
  1798.     r=ground(m[i1][i]);f=subri(m[i1][i],r);
  1799.     p3=gsub(b[i1],gmul(r,b[i]));
  1800.     b[i1]=b[i];b[i]=p3;
  1801.     for(j=1;j<i;j++)
  1802.     {
  1803.       if(!qzer[j])
  1804.       {
  1805.         p3=mpsub(m[i1][j],mulir(r,m[i][j]));
  1806.         affrr(m[i][j],m[i1][j]);mpaff(p3,m[i][j]);
  1807.       }
  1808.     }
  1809.     c1=addrr(bn[i1],mulrr(mulrr(f,f),bn[i]));
  1810.     fl=quazero(c1);
  1811.     if(!fl)
  1812.     {
  1813.       c2=divrr(mulrr(bn[i],f),c1);affrr(c2,m[i1][i]);
  1814.       c3=divrr(bn[i1],c1);mulrrz(c3,bn[i],bn[i1]);
  1815.       affrr(c1,bn[i]);qzer[i1]=quazero(bn[i1]);qzer[i]=0;
  1816.       for(j=i+2;j<=n;j++)
  1817.       {
  1818.         p3=addrr(mulrr(m[j][i1],c3),mulrr(m[j][i],c2));
  1819.         subrrz(m[j][i],mulrr(f,m[j][i1]),m[j][i1]);
  1820.         affrr(p3,m[j][i]);
  1821.       }
  1822.     }
  1823.     else
  1824.     {
  1825.       qzer[i1]=qzer[i];affrr(bn[i],bn[i1]);
  1826.       affrr(c1,bn[i]);qzer[i]=1;
  1827.       for(j=i+2;j<=n;j++) affrr(m[j][i],m[j][i1]);
  1828.     }
  1829.   }
  1830.   p3=cgetg(lx,18);
  1831.   for(i=1;i<=n;i++) p3[i]=(i==n) ? un : zero;
  1832.   p3[n]=un;p5=cgetg(lx,19);
  1833.   for(i=1;i<=n;i++) p5[i]=lcopy(b[i]);
  1834.   p4=gauss(gtrans(p5),p3);tetpil=avma;
  1835.   y=gerepile(av,tetpil,gtrans(p4));
  1836.   return y;
  1837. }
  1838.  
  1839. GEN    algdep(x,n,prec)
  1840.      
  1841.      GEN  x;
  1842.      long n,prec;
  1843.      
  1844. {
  1845.   long   tx=typ(x),av,tetpil,i,j,k;
  1846.   GEN    y,p1;
  1847.   
  1848.   if(tx>=10) err(algder1);
  1849.   if(tx==9) {y=gcopy(x[1]);setvarn(y,0);return y;}
  1850.   if(gcmp0(x)) return gzero;
  1851.   av=avma;p1=cgetg(n+2,18);p1[1]=un;
  1852.   for(i=2;i<=n+1;i++) p1[i]=lmul(p1[i-1],x);
  1853.   p1=lindep(p1,prec);
  1854.   tetpil=avma;y=cgetg(n+3,10);
  1855.   setsigne(y,1);setvarn(y,0);j=0;
  1856.   k=1;while(gcmp0(p1[k])) k++;
  1857.   for(i=0;i<=n+1-k;i++)
  1858.   {
  1859.     y[i+2]=lcopy(p1[k+i]);
  1860.     if(!gcmp0(p1[k+i])) j=i;
  1861.   }
  1862.   setlgef(y,j+3);
  1863.   if(gsigne(y[j+2])>0) return gerepile(av,tetpil,y);
  1864.   else {tetpil=avma;return gerepile(av,tetpil,gneg(y));}
  1865. }
  1866.  
  1867. GEN    algdep2(x,n,bit,prec)
  1868.      
  1869.      GEN  x;
  1870.      long n,bit,prec;
  1871.      
  1872. {
  1873.   long   tx=typ(x),av,tetpil,i,j,k;
  1874.   GEN    y,p1;
  1875.   
  1876.   if(tx>=10) err(algder1);
  1877.   if(tx==9) {y=gcopy(x[1]);setvarn(y,0);return y;}
  1878.   if(gcmp0(x)) return gzero;
  1879.   av=avma;p1=cgetg(n+2,18);p1[1]=un;
  1880.   for(i=2;i<=n+1;i++) p1[i]=lmul(p1[i-1],x);
  1881.   p1=lindep2(p1,bit,prec);
  1882.   tetpil=avma;y=cgetg(n+3,10);
  1883.   setsigne(y,1);setvarn(y,0);j=0;
  1884.   k=1;while(gcmp0(p1[k])) k++;
  1885.   for(i=0;i<=n+1-k;i++)
  1886.   {
  1887.     y[i+2]=lcopy(p1[k+i]);
  1888.     if(!gcmp0(p1[k+i])) j=i;
  1889.   }
  1890.   setlgef(y,j+3);
  1891.   if(gsigne(y[j+2])>0) return gerepile(av,tetpil,y);
  1892.   else {tetpil=avma;return gerepile(av,tetpil,gneg(y));}
  1893. }
  1894.  
  1895. /********************************************************************/
  1896. /********************************************************************/
  1897. /**                                                                **/
  1898. /**                   CHANGEMENTS DE VARIABLES                     **/
  1899. /**                                                                **/
  1900. /********************************************************************/
  1901. /********************************************************************/
  1902.  
  1903. GEN  changevar(x,y)
  1904.      
  1905.      GEN  x,y;
  1906.      
  1907. /* Substitution globale des composantes du vecteur y aux variables de x */
  1908.   
  1909. {
  1910.   long tx=typ(x),ty=typ(y),lx=lg(x),ly=lg(y),vx,vy,i,av,tetpil;
  1911.   GEN  p1,p2,p3,z;
  1912.   
  1913.   if((ty<17) || (ty>18)) err(changer1);
  1914.   if(tx<9) return gcopy(x);
  1915.   if(tx==9)
  1916.     {
  1917.       av=avma;p1=changevar(x[1],y);p2=changevar(x[2],y);
  1918.       if(isonstack(x[1])) 
  1919.     {tetpil=avma;return gerepile(av,tetpil,gmodulcp(p1,p2));}
  1920.       else {tetpil=avma;return gerepile(av,tetpil,gmodulo(p1,p2));}
  1921.     }
  1922.   if((tx==13)||(tx==14))
  1923.     {
  1924.       av=avma;p1=changevar(x[1],y);p2=changevar(x[2],y);
  1925.       tetpil=avma;return gerepile(av,tetpil,gdiv(p1,p2));
  1926.     }
  1927.   if(tx>11)
  1928.   {
  1929.     z=cgetg(lx,tx);
  1930.     for(i=1;i<lx;i++) z[i]=lchangevar(x[i],y);
  1931.     return z;
  1932.   }
  1933.   vx=varn(x)+1;if(vx>=ly) return gcopy(x);
  1934.   if(!signe(x))
  1935.   {
  1936.     vy=gvar(y[vx]); if(vy>255) err(changer1);
  1937.     z=gcopy(x);
  1938.     setvarn(z,vy);
  1939.   }
  1940.   else
  1941.   {
  1942.     av=avma;p1=(GEN)y[vx];
  1943.     if(tx==10)
  1944.     {
  1945.       lx=lgef(x);
  1946.       p2=changevar(x[lx-1],y);
  1947.       for(i=lx-2;i>=2;i--)
  1948.       {
  1949.         p2=gmul(p2,p1);p3=changevar(x[i],y);
  1950.         tetpil=avma;p2=gadd(p2,p3);
  1951.       }
  1952.       if(lx>3) z=gerepile(av,tetpil,p2);
  1953.       else z=p2;
  1954.     }
  1955.     else
  1956.     {
  1957.       p2=changevar(x[lx-1],y);
  1958.       for(i=lx-2;i>=2;i--)
  1959.       {
  1960.         p2=gmul(p2,p1);p3=changevar(x[i],y);
  1961.         p2=gadd(p2,p3);
  1962.       }
  1963.       p3=ggrando(p1,lx-2);tetpil=avma;
  1964.       z=gadd(p2,p3);
  1965.       if(valp(x))
  1966.       {
  1967.         p2=gpuigs(p1,valp(x));tetpil=avma;
  1968.         z=gerepile(av,tetpil,gmul(p2,z));
  1969.       }
  1970.       else z=gerepile(av,tetpil,z);
  1971.     }
  1972.   }
  1973.   return z;
  1974. }
  1975.  
  1976. GEN reorder(x)
  1977.   GEN x;
  1978.      
  1979. {
  1980.   long tx=typ(x), lx=lg(x), t1[MAXVAR], i, j, n;
  1981.   
  1982.   if((tx < 17) || (tx > 18)) err(reorder1);
  1983.   for(i=0; i<nvar; i++) t1[i] = 0;
  1984.   for(n=1; n<lx; n++)
  1985.   {
  1986.     i = gvar(x[n]);
  1987.     if (i >= nvar) err(reorder2);
  1988.     if (t1[i]) err(reorder3);
  1989.     t1[i] = 1;
  1990.   }
  1991.   for(i=n=0; i<nvar; i++)
  1992.     if (t1[i])
  1993.     {
  1994.       j = gvar(x[++n]);
  1995.       polvar[i+1]=lpolx[j];
  1996.       ordvar[j]=i;
  1997.     }
  1998.   varchanged=0;
  1999.   for(i=0;i<nvar;i++) if(ordvar[i]!=i) {varchanged=1;break;}
  2000.   return polvar;
  2001. }
  2002.  
  2003. /********************************************************************/
  2004. /********************************************************************/
  2005. /**                                                                **/
  2006. /**                       TRI PAR HEAPSORT                         **/
  2007. /**                                                                **/
  2008. /********************************************************************/
  2009. /********************************************************************/
  2010.  
  2011. GEN sort(x)
  2012.  
  2013.      GEN x;
  2014.  
  2015. {
  2016.   long av,tetpil,n,i,j,indxt,q,ir,l,tx=typ(x),lx=lg(x);
  2017.   GEN y,indx;
  2018.  
  2019.   if((tx<17)||(tx>18)) err(sorter1);
  2020.   av=avma;n=lx-1;if(n<=1) return gcopy(x);
  2021.   indx=cgeti(lx);for(j=1;j<=n;j++) indx[j]=j;
  2022.   l=(n>>1)+1;ir=n;
  2023.   for(;;)
  2024.   {
  2025.     if(l>1) q=x[(indxt=indx[--l])];
  2026.     else
  2027.     {
  2028.       q=x[(indxt=indx[ir])];indx[ir]=indx[1];
  2029.       if(--ir ==1)
  2030.       {
  2031.         indx[1]=indxt;tetpil=avma;y=cgetg(lx,tx);
  2032.         for(i=1;i<=n;i++) y[i]=lcopy(x[indx[i]]);
  2033.         return gerepile(av,tetpil,y);
  2034.       }
  2035.     }
  2036.     i=l;j=l<<1;
  2037.     while(j<=ir)
  2038.     {
  2039.       if(j<ir && gcmp(x[indx[j]],x[indx[j+1]])<0) j++;
  2040.       if(gcmp(q,x[indx[j]])<0) {indx[i]=indx[j];j+=(i=j);}
  2041.       else j=ir+1;
  2042.     }
  2043.     indx[i]=indxt;
  2044.   }
  2045. }
  2046.  
  2047. GEN lexsort(x)
  2048.  
  2049.      GEN x;
  2050.  
  2051. {
  2052.   long av,tetpil,n,i,j,indxt,q,ir,l,tx=typ(x),lx=lg(x);
  2053.   GEN y,indx;
  2054.  
  2055.   if((tx<17)||(tx>18)) err(sorter1);
  2056.   av=avma;n=lx-1;if(n<=1) return gcopy(x);
  2057.   indx=cgeti(lx);for(j=1;j<=n;j++) indx[j]=j;
  2058.   l=(n>>1)+1;ir=n;
  2059.   for(;;)
  2060.   {
  2061.     if(l>1) q=x[(indxt=indx[--l])];
  2062.     else
  2063.     {
  2064.       q=x[(indxt=indx[ir])];indx[ir]=indx[1];
  2065.       if(--ir ==1)
  2066.       {
  2067.         indx[1]=indxt;tetpil=avma;y=cgetg(lx,tx);
  2068.         for(i=1;i<=n;i++) y[i]=lcopy(x[indx[i]]);
  2069.         return gerepile(av,tetpil,y);
  2070.       }
  2071.     }
  2072.     i=l;j=l<<1;
  2073.     while(j<=ir)
  2074.     {
  2075.       if(j<ir && lexcmp(x[indx[j]],x[indx[j+1]])<0) j++;
  2076.       if(lexcmp(q,x[indx[j]])<0) {indx[i]=indx[j];j+=(i=j);}
  2077.       else j=ir+1;
  2078.     }
  2079.     indx[i]=indxt;
  2080.   }
  2081. }
  2082.  
  2083.  
  2084. GEN vecsort(x,k)
  2085.      GEN x;
  2086.      long k;
  2087. {
  2088.   long av,tetpil,n,i,j,indxt,q,ir,l,t,tx=typ(x),lx=lg(x);
  2089.   GEN y,indx;
  2090.  
  2091.   if(tx<17) err(sorter1);
  2092.   av=avma;n=lx-1;if(n<=1) return gcopy(x);
  2093.   indx=cgeti(lx);
  2094.   for(j=1;j<=n;j++)
  2095.     {
  2096.       indx[j]=j;t=typ(x[j]);
  2097.       if((t<17)||(t>18)) err(sorter1);
  2098.     }
  2099.   l=(n>>1)+1;ir=n;
  2100.   for(;;)
  2101.   {
  2102.     if(l>1) q=x[(indxt=indx[--l])];
  2103.     else
  2104.     {
  2105.       q=x[(indxt=indx[ir])];indx[ir]=indx[1];
  2106.       if(--ir ==1)
  2107.       {
  2108.         indx[1]=indxt;tetpil=avma;y=cgetg(lx,tx);
  2109.         for(i=1;i<=n;i++) y[i]=lcopy(x[indx[i]]);
  2110.         return gerepile(av,tetpil,y);
  2111.       }
  2112.     }
  2113.     i=l;j=l<<1;
  2114.     while(j<=ir)
  2115.     {
  2116.       if(j<ir && gcmp(((GEN)x[indx[j]])[k],((GEN)x[indx[j+1]])[k])<0) j++;
  2117.       if(gcmp(((GEN)q)[k],((GEN)x[indx[j]])[k])<0) {indx[i]=indx[j];j+=(i=j);}
  2118.       else j=ir+1;
  2119.     }
  2120.     indx[i]=indxt;
  2121.   }
  2122. }
  2123.  
  2124. GEN indexsort(x)
  2125.  
  2126.      GEN x;
  2127.  
  2128. {
  2129.   long av,tetpil,n,i,j,indxt,q,ir,l,tx=typ(x),lx=lg(x);
  2130.   GEN y,indx;
  2131.  
  2132.   if((tx<17)||(tx>18)) err(sorter1);
  2133.   av=avma;n=lx-1;if(n<=1) {y=cgetg(lx,tx);y[1]=un;return y;}
  2134.   indx=cgeti(lx);for(j=1;j<=n;j++) indx[j]=j;
  2135.   l=(n>>1)+1;ir=n;
  2136.   for(;;)
  2137.   {
  2138.     if(l>1) q=x[(indxt=indx[--l])];
  2139.     else
  2140.     {
  2141.       q=x[(indxt=indx[ir])];indx[ir]=indx[1];
  2142.       if(--ir ==1)
  2143.       {
  2144.         indx[1]=indxt;tetpil=avma;y=cgetg(lx,tx);
  2145.         for(i=1;i<=n;i++) y[i]=lstoi(indx[i]);
  2146.         return gerepile(av,tetpil,y);
  2147.       }
  2148.     }
  2149.     i=l;j=l<<1;
  2150.     while(j<=ir)
  2151.     {
  2152.       if(j<ir && gcmp(x[indx[j]],x[indx[j+1]])<0) j++;
  2153.       if(gcmp(q,x[indx[j]])<0) {indx[i]=indx[j];j+=(i=j);}
  2154.       else j=ir+1;
  2155.     }
  2156.     indx[i]=indxt;
  2157.   }
  2158. }
  2159.  
  2160. /********************************************************************/
  2161. /********************************************************************/
  2162. /**                                                                **/
  2163. /**                  INTERPOLATION POLYNOMIALE                     **/
  2164. /**                                                                **/
  2165. /********************************************************************/
  2166. /********************************************************************/
  2167.  
  2168. GEN  polint(xa,ya,x,dy)
  2169.  
  2170.      GEN xa,ya,x,*dy;
  2171.  
  2172. {
  2173.   long av,av1,tetpil,dec,i,m,ns=1,tx=typ(xa),ty=typ(ya),n,lx=lg(xa),ly=lg(ya);
  2174.   GEN den,dif,dift,ho,hp,w,y,c,d;
  2175.  
  2176.   if((tx<17)||(tx>18)||(ty<17)||(ty>18)||(lx!=ly)) err(polinter1);
  2177.   n=lx-1;if(n<=1) {y=gcopy(ya[1]);*dy=gcopy(y);return y;}
  2178.   av=avma;c=cgetg(ly,ty);d=cgetg(ly,ty);
  2179.   dif=gabs(gsub(x,xa[1]));
  2180.   for(i=1;i<=n;i++)
  2181.     {
  2182.       if(gcmp((dift=gabs(gsub(x,xa[i]))),dif)<0) {ns=i;dif=dift;}
  2183.       c[i]=ya[i];d[i]=ya[i];
  2184.     }
  2185.   y=(GEN)ya[ns--];
  2186.   for(m=1;m<n;m++)
  2187.     {
  2188.       for(i=1;i<=n-m;i++)
  2189.         {
  2190.           ho=gsub(xa[i],x);hp=gsub(xa[i+m],x);w=gsub(c[i+1],d[i]);
  2191.           if(gcmp0(den=gsub(ho,hp))) err(polinter2);
  2192.           den=gdiv(w,den);d[i]=lmul(hp,den);c[i]=lmul(ho,den);
  2193.         }
  2194.       *dy=(2*ns<(n-m))?(GEN)c[ns+1]:(GEN)d[ns--];
  2195.       tetpil=avma;y=gadd(y,*dy);
  2196.     }
  2197.   *dy=gcopy(*dy);av1=avma;dec=lpile(av,tetpil,0)>>2;
  2198.   if(adecaler(y,tetpil,av1)) y+=dec;
  2199.   if(adecaler(*dy,tetpil,av1)) (*dy)+=dec;
  2200.   return y;
  2201. }
  2202.  
  2203. /********************************************************************/
  2204. /********************************************************************/
  2205. /**                                                                **/
  2206. /**                    POLRED ET COMPAGNIE                         **/
  2207. /**                                                                **/
  2208. /********************************************************************/
  2209. /********************************************************************/
  2210.  
  2211. GEN polsym(x,n)
  2212.      GEN x;
  2213.      long n;
  2214. {
  2215.   long av1,av2,tx=typ(x),dx=lgef(x)-3,i,k;
  2216.   GEN s,y;
  2217.  
  2218.   if((tx!=10)||(!signe(x))) err(poltyper);
  2219.   y=cgetg(n+2,18);y[1]=lstoi(dx);
  2220.   if(n<0) err(impl,"polsym of a negative n");
  2221.   for(k=1;k<=n;k++)
  2222.     {
  2223.       av1=avma;s=(dx>=k) ? gmulsg(k,x[dx+2-k]) : gzero;
  2224.       for(i=1;(i<k)&&(i<=dx);i++)
  2225.     s=gadd(s,gmul(y[k-i+1],x[dx+2-i]));
  2226.       if(!gcmp1(x[dx+2])) s=gdiv(s,x[dx+2]);
  2227.       av2=avma;y[k+1]=lpile(av1,av2,gneg(s));
  2228.     }
  2229.   return y;
  2230. }
  2231.  
  2232. GEN allpolred(x,pta,code,prec)
  2233.      GEN x,*pta;
  2234.      long code,prec;
  2235. {
  2236.   GEN y,p1,p2,p3,p4,p5,p6,p7,ptrace,a;
  2237.   long tx=typ(x),n=lgef(x)-3,i,j,k,dec,av=avma,av1,tet1,tet2,tetpil,v;
  2238.  
  2239.   if(tx!=10) err(poltyper);
  2240.   if(!signe(x)) return gcopy(x); 
  2241.   if(!gcmp1(x[n+2])) err(poltyper);
  2242.   p4=allbase(x,code,&p2);
  2243.   if(sturm(x)<n)
  2244.     {
  2245.       p2=roots(x,prec);p3=cgetg(n+1,19);
  2246.       for(i=1;i<=n;i++)
  2247.     {
  2248.       p1=cgetg(n+1,18);p3[i]=(long)p1;
  2249.       for(j=1;j<=n;j++)
  2250.         p1[j]=lsubst(p4[i],varn(p4[i]),p2[j]);
  2251.     }
  2252.       p2=greal(gmul(gconj(gtrans(p3)),p3));
  2253.     }
  2254.   else
  2255.     {
  2256.       ptrace=cgetg(n+1,17);ptrace[1]=lstoi(n);
  2257.       for(k=1;k<n;k++) 
  2258.     {
  2259.       p3=gmulsg(k,x[n-k+2]);
  2260.       for(i=1;i<k;i++) p3=gadd(p3,gmul(x[n-i+2],ptrace[k-i+1]));
  2261.       ptrace[k+1]=lneg(p3);
  2262.     }
  2263.       p2=cgetg(n+1,19);
  2264.       for(i=1;i<=n;i++)
  2265.     {
  2266.       p1=cgetg(n+1,18);p2[i]=(long)p1;
  2267.       for(j=1;j<i;j++) p1[j]=lcopy(coeff(p2,i,j));
  2268.       for(j=i;j<=n;j++)
  2269.         {
  2270.           p5=gres(gmul(p4[i],p4[j]),x);p6=gzero;
  2271.           for(k=0;k<=lgef(p5)-3;k++) p6=gadd(p6,gmul(p5[k+2],ptrace[k+1]));
  2272.           p1[j]=(long)p6;
  2273.         }
  2274.     }
  2275.     }
  2276.   p1=lllgram(p2,prec);v=varn(x);tet1=avma;
  2277.   a=cgetg(n+1,18);for(i=1;i<=n;i++) a[i]=lmul(p4,p1[i]);
  2278.   tetpil=avma;y=cgetg(n+1,18);
  2279.   for(i=1;i<=n;i++)
  2280.     {
  2281.       av1=avma;p3=gmodulcp(a[i],x);p7=content(p3[2]);
  2282.       if(gcmp1(p7)) p3=caract(p3,v);
  2283.       else
  2284.     {
  2285.       p3=caract(gdiv(p3,p7),v);
  2286.       p3=gmul(gpuigs(p7,lgef(p3)-3),gsubst(p3,v,gdiv(polx[v],p7)));
  2287.     }
  2288.       p5=ggcd(deriv(p3,v),p3);
  2289.       p6=(GEN)p5[lgef(p5)-1];if(!gcmp1(p6)) p5=gdiv(p5,p6);
  2290.       tet2=avma;p3=gerepile(av1,tet2,gdiv(p3,p5));
  2291.       y[i]=(long)p3;j=lgef(p3)-2;
  2292.       while((j>=2)&&(!signe(p3[j]))) j-=2;
  2293.       if((j>=2)&&(signe(p3[j])>0))
  2294.     {
  2295.       for(;j>=2;j-=2) setsigne(p3[j],-signe(p3[j]));
  2296.       gnegz(a[i],a[i]);
  2297.     }
  2298.     }
  2299.   if(pta) {dec=lpile(av,tet1,0)>>2;*pta=a+dec;return y+dec;}
  2300.   else return gerepile(av,tetpil,y);
  2301. }
  2302.  
  2303. GEN polred(x,prec)
  2304.      GEN x;
  2305.      long prec;
  2306. {
  2307.   return allpolred(x,(GEN)0,0,prec);
  2308. }
  2309.  
  2310. GEN smallpolred(x,prec)
  2311.      GEN x;
  2312.      long prec;
  2313. {
  2314.   return allpolred(x,(GEN)0,1,prec);
  2315. }
  2316.  
  2317. GEN factoredpolred(x,p,prec)
  2318.      GEN x,p;
  2319.      long prec;
  2320. {
  2321.   return allpolred(x,(GEN)0,(long)p,prec);
  2322. }
  2323.  
  2324. GEN polred2(x,prec)
  2325.      GEN x;
  2326.      long prec;
  2327. {
  2328.   GEN y;
  2329.  
  2330.   y=cgetg(3,19);
  2331.   y[2]=(long)allpolred(x,y+1,0,prec);
  2332.   return y;
  2333. }
  2334.  
  2335. GEN smallpolred2(x,prec)
  2336.      GEN x;
  2337.      long prec;
  2338. {
  2339.   GEN y;
  2340.  
  2341.   y=cgetg(3,19);
  2342.   y[2]=(long)allpolred(x,y+1,1,prec);
  2343.   return y;
  2344. }
  2345.  
  2346. GEN factoredpolred2(x,p,prec)
  2347.      GEN x,p;
  2348.      long prec;
  2349. {
  2350.   GEN y;
  2351.  
  2352.   y=cgetg(3,19);
  2353.   y[2]=(long)allpolred(x,y+1,(long)p,prec);
  2354.   return y;
  2355. }
  2356.  
  2357. GEN ordred(x,prec)
  2358.      GEN x;
  2359.      long prec;
  2360. {
  2361.   GEN y,p1,p2,p3,p4,p5,p6,p7;
  2362.   long tx=typ(x),n=lgef(x)-3,i,j,av=avma,v,av1,tet1,tet2;
  2363.  
  2364.   if(tx!=10) err(poltyper);
  2365.   if(!signe(x)) return gcopy(x); 
  2366.   if(!gcmp1(x[n+2])) err(poltyper);
  2367.   p2=roots(x,prec);p3=cgetg(n+1,19);
  2368.   p4=cgetg(n+1,17);for(i=1;i<=n;i++) p4[i]=lpuigs(polx[varn(x)],i-1);
  2369.   for(i=1;i<=n;i++)
  2370.     {
  2371.       p1=cgetg(n+1,18);p3[i]=(long)p1;
  2372.       for(j=1;j<=n;j++)
  2373.     p1[j]=lpuigs(p2[j],i-1);
  2374.     }
  2375.   p2=greal(gmul(gconj(gtrans(p3)),p3));
  2376.   p1=lllgram(p2,prec);v=varn(x);tet1=avma;y=cgetg(n+1,18);
  2377.   for(i=1;i<=n;i++)
  2378.     {
  2379.       av1=avma;p3=gmodulcp(gmul(p4,p1[i]),x);p7=content(p3[2]);
  2380.       if(gcmp1(p7)) p3=caract(p3,v);
  2381.       else
  2382.     {
  2383.       p3=caract(gdiv(p3,p7),v);
  2384.       p3=gmul(gpuigs(p7,lgef(p3)-3),gsubst(p3,v,gdiv(polx[v],p7)));
  2385.     }
  2386.       p5=ggcd(deriv(p3,v),p3);
  2387.       p6=(GEN)p5[lgef(p5)-1];if(!gcmp1(p6)) p5=gdiv(p5,p6);
  2388.       tet2=avma;p3=gerepile(av1,tet2,gdiv(p3,p5));
  2389.       y[i]=(long)p3;j=lgef(p3)-2;
  2390.       while((j>=2)&&(!signe(p3[j]))) j-=2;
  2391.       if((j>=2)&&(signe(p3[j])>0))
  2392.     {for(;j>=2;j-=2) setsigne(p3[j],-signe(p3[j]));}
  2393.     }
  2394.   return gerepile(av,tet1,y);
  2395. }
  2396.  
  2397. /*===================================================
  2398.   Nombre de vecteurs minimaux du reseau defini par
  2399.               la matrice de GRAM  a
  2400. ====================================================*/
  2401.  
  2402. GEN minim1(a)
  2403.      GEN a;
  2404.  
  2405. {
  2406.   GEN u,r,unr;
  2407.   long n1=lg(a),n,av,nonnul,i,j,k,s,borne,norme,*x;
  2408.   double eps=0.000001,p,b,c;
  2409.   double **q,*v,*y,*z;
  2410.  
  2411.   n=n1-1;
  2412.   x=(long*)malloc(n1*sizeof(long));
  2413.   y=(double*)malloc(n1*sizeof(double));
  2414.   z=(double*)malloc(n1*sizeof(double));
  2415.   v=(double*)malloc(n1*sizeof(double));
  2416.   q=(double**)malloc(4*n1);
  2417.   for(j=1;j<=n;j++) q[j]=(double*)malloc(n1*sizeof(double));
  2418.  
  2419.   av=avma;
  2420.   affsr(1,unr=cgetr(6));
  2421.   u=lllgram(a,6);
  2422.   a=gmul(a,unr);a=gmul(gtrans(u),gmul(a,u));
  2423.   r=sqred1(a);
  2424.   for(j=1;j<=n;j++)
  2425.     {
  2426.       v[j]=rtodbl(coeff(r,j,j));
  2427.       for(i=1;i<j;i++)
  2428.     q[i][j]=rtodbl(coeff(r,i,j));
  2429.     }
  2430.   b=rtodbl(coeff(a,1,1));
  2431.   for(i=2;i<=n;i++)
  2432.     if((c=rtodbl(coeff(a,i,i)))<b) b=c;
  2433.   avma=av;
  2434.   borne=b+eps;
  2435.   s=0;
  2436.   do
  2437.     {
  2438.       k=n;
  2439.       y[n]=z[n]=0;
  2440.       x[n]=sqrt(borne/v[n]+eps);
  2441.       do
  2442.     {
  2443.       do
  2444.         {
  2445.           if(k>1)
  2446.         {
  2447.           k--;
  2448.           z[k]=0;
  2449.           for(j=k+1;j<=n;j++) z[k]=z[k]+q[k][j]*x[j];
  2450.           p=x[k+1]+z[k+1];
  2451.           y[k]=y[k+1]+p*p*v[k+1];
  2452.           x[k]=floor(sqrt((borne-y[k]+eps)/v[k])-z[k]);
  2453.         }
  2454.           while(v[k]*(x[k]+z[k])*(x[k]+z[k])>borne-y[k]+eps)
  2455.         {k++;x[k]--;}
  2456.         }
  2457.       while(k>1);
  2458.       if(nonnul=(x[1]||(y[1]>eps)))
  2459.         {
  2460.           norme=y[1]+v[1]*(x[1]+z[1])*(x[1]+z[1])+eps;
  2461.           if (norme==borne)    { s++;x[k]--;}
  2462.           else
  2463.         { s=0; borne=norme;}
  2464.         }
  2465.     }
  2466.       while(nonnul && s);
  2467.     }
  2468.   while(nonnul);
  2469.   free(x);free(y);free(z);
  2470.   for(j=1;j<=n;j++) free(q[j]);free(q);
  2471.   u=cgetg(3,17);
  2472.   u[1]=lstoi(s<<1);u[2]=lstoi(borne);return u;
  2473. }
  2474.  
  2475. GEN minim(a,borne,stockmax)
  2476.      GEN a;
  2477.      long borne,stockmax;
  2478.  
  2479. {
  2480.   GEN u,r,unr,S,S1;
  2481.   long n1=lg(a),n,av,av1,nonnul,i,j,k,s,norme,flg,*x;
  2482.   double p,b,c;
  2483.   double **q,*v,*y,*z;
  2484.   double eps=0.000001;
  2485.  
  2486.   n=n1-1;
  2487.   x=(long*)malloc(n1*sizeof(long));
  2488.   y=(double*)malloc(n1*sizeof(double));
  2489.   z=(double*)malloc(n1*sizeof(double));
  2490.   v=(double*)malloc(n1*sizeof(double));
  2491.   q=(double**)malloc(4*n1);
  2492.   for(j=1;j<=n;j++) q[j]=(double*)malloc(n1*sizeof(double));
  2493.   av=avma;
  2494.   affsr(1,unr=cgetr(6));
  2495.   u=lllgram(a,6);
  2496.   a=gmul(a,unr);a=gmul(gtrans(u),gmul(a,u));
  2497.   r=sqred1(a);
  2498.   for(j=1;j<=n;j++)
  2499.     {
  2500.       v[j]=rtodbl(coeff(r,j,j));
  2501.       for(i=1;i<j;i++)
  2502.     q[i][j]=rtodbl(coeff(r,i,j));
  2503.     }
  2504.    if(!(flg=borne))   /* flg=0 <==> chercher les vec min */
  2505.      {
  2506.        b=rtodbl(coeff(a,1,1));
  2507.        for(i=2;i<=n;i++)
  2508.      if((c=rtodbl(coeff(a,i,i)))<b) b=c;
  2509.        borne=b+eps;
  2510.      }
  2511.   if(stockmax) S=cgetg(stockmax+1,19);
  2512.   s=0;k=n;
  2513.   y[n]=z[n]=0;
  2514.   x[n]=sqrt(borne/v[n]+eps);
  2515.   do
  2516.     {
  2517.       do
  2518.     {
  2519.       if(k>1)
  2520.         {
  2521.           k--;
  2522.           z[k]=0;
  2523.           for(j=k+1;j<=n;j++) z[k]=z[k]+q[k][j]*x[j];
  2524.           p=x[k+1]+z[k+1];
  2525.           y[k]=y[k+1]+p*p*v[k+1];
  2526.           x[k]=floor(sqrt((borne-y[k]+eps)/v[k])-z[k]);
  2527.         }
  2528.       while(v[k]*(x[k]+z[k])*(x[k]+z[k])>borne-y[k]+eps)
  2529.         {k++;x[k]--;}
  2530.     }
  2531.       while(k>1);
  2532.       if(nonnul=(x[1]||(y[1]>eps)))
  2533.     {
  2534.       norme=y[1]+v[1]*(x[1]+z[1])*(x[1]+z[1])+eps;
  2535.       if(!flg&&(norme<borne))
  2536.         { s=0; borne=norme;};
  2537.       s++;
  2538.       if(s<=stockmax)
  2539.         {
  2540.           S[s]=lgetg(n+1,18);
  2541.           for(i=1;i<=n;i++) coeff(S,i,s)=lstoi(x[i]);
  2542.         }
  2543.       x[k]--;
  2544.     }
  2545.     }
  2546.   while(nonnul);
  2547.   free(x);free(y);free(z);free(v);
  2548.   for(j=1;j<=n;j++) free(q[j]);free(q);
  2549.   if(stockmax) 
  2550.     {
  2551.       av1=avma;
  2552.       k=(s<stockmax)? s:stockmax;
  2553.       S1=cgetg(k+1,19);for(j=1;j<=k;j++) S1[j]=lmul(u,S[j]);
  2554.       S=gerepile(av,av1,S1);
  2555.     }
  2556.   else {avma=av;S=cgetg(1,19);}
  2557.   u=cgetg(4,17);
  2558.   u[1]=lstoi(s<<1);
  2559.   u[2]=lstoi(norme);
  2560.   u[3]=(long)S;
  2561.   return u;
  2562. }
  2563.  
  2564. GEN polymodrecip(x)
  2565.      GEN x;
  2566. {
  2567.   long v,i,j,n,av,tetpil,lx;
  2568.   GEN p1,p2,p3,p,phi,y,col;
  2569.  
  2570.   if(typ(x)!=9) err(polymoder1);
  2571.   p=(GEN)x[1];phi=(GEN)x[2];if(gcmp0(phi)) err(polymoder2);
  2572.   v=varn(p);n=lgef(p)-3;if(n<=0) return gcopy(x);
  2573.   if(n==1)
  2574.     {
  2575.       y=cgetg(3,9);p1=cgetg(4,10);y[1]=(long)p1;
  2576.       p1[1]=p[1];p1[2]=(typ(phi)==10) ? lneg(phi[2]) : lneg(phi);
  2577.       p1[3]=un;p1=cgetg(3,10);y[2]=(long)p1;
  2578.       if(gcmp0(p[2])) p1[1]=2;
  2579.       else
  2580.     {
  2581.       p1[1]=p[1]-1;av=avma;p2=gdiv(p[2],p[3]);
  2582.       tetpil=avma;p1[2]=lpile(av,tetpil,gneg(p2));
  2583.     }
  2584.       setvarn(p1,v);
  2585.       return y;
  2586.     }
  2587.   else
  2588.     {
  2589.       av=avma;y=cgetg(n+1,19);p1=cgetg(n+1,18);
  2590.       y[1]=(long)p1;p1[1]=un;for(i=2;i<=n;i++) p1[i]=zero;
  2591.       p2=phi;
  2592.       for(j=2;j<=n;j++)
  2593.     {
  2594.       lx=lgef(p2);p1=cgetg(n+1,18);y[j]=(long)p1;
  2595.       for(i=1;i<=lx-2;i++) p1[i]=p2[i+1];
  2596.       for(i=lx-1;i<=n;i++) p1[i]=zero;
  2597.       if(j<n) p2=gmod(gmul(p2,phi),p);
  2598.     }
  2599.       col=cgetg(n+1,18);col[1]=zero;col[2]=un;
  2600.       for(i=3;i<=n;i++) col[i]=zero;
  2601.       p1=gauss(y,col);p2=gtopolyrev(p1,v);
  2602.       p3=caract(x,v);
  2603.       tetpil=avma;return gerepile(av,tetpil,gmodulcp(p2,p3));
  2604.     }
  2605. }
  2606.  
  2607. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  2608. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  2609. /*~                                    ~*/
  2610. /*~                     RANDOM                    ~*/
  2611. /*~                                    ~*/
  2612. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  2613. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  2614.  
  2615. GEN genrand()
  2616. {
  2617.   return stoi(rand());
  2618. }
  2619.  
  2620.  
  2621. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  2622. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  2623. /*~                                    ~*/
  2624. /*~                 PERMUTATIONS                ~*/
  2625. /*~                                    ~*/
  2626. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  2627. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  2628.  
  2629. GEN permute(n,x)
  2630.      long n;
  2631.      GEN x;
  2632. {
  2633.   long av=avma,tetpil,i,a,r;
  2634.   GEN v,w,y;
  2635.  
  2636.   v=cgetg(2,17);v[1]=1;
  2637.   for(r=2;r<=n;r++)
  2638.     {
  2639.       x=dvmdis(x,r,&w);a=itos(w);
  2640.       w=cgetg(r+1,17);for(i=1;i<=a;i++) w[i]=v[i];
  2641.       w[a+1]=r;for(i=a+2;i<=r;i++) w[i]=v[i-1];
  2642.       v=w;
  2643.     }
  2644.   tetpil=av;y=cgetg(n+1,17);
  2645.   for(i=1;i<=n;i++) y[i]=lstoi(v[i]);
  2646.   return gerepile(av,tetpil,y);
  2647. }
  2648.