===================================================================
RCS file: /home/cvs/OpenXM_contrib/pari-2.2/src/basemath/Attic/bibli2.c,v
retrieving revision 1.1
retrieving revision 1.2
diff -u -p -r1.1 -r1.2
--- OpenXM_contrib/pari-2.2/src/basemath/Attic/bibli2.c	2001/10/02 11:17:02	1.1
+++ OpenXM_contrib/pari-2.2/src/basemath/Attic/bibli2.c	2002/09/11 07:26:50	1.2
@@ -1,4 +1,4 @@
-/* $Id: bibli2.c,v 1.1 2001/10/02 11:17:02 noro Exp $
+/* $Id: bibli2.c,v 1.2 2002/09/11 07:26:50 noro Exp $
 
 Copyright (C) 2000  The PARI group.
 
@@ -20,6 +20,7 @@ Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 /**                                                                **/
 /********************************************************************/
 #include "pari.h"
+extern GEN _ei(long n, long i);
 
 /********************************************************************/
 /**                                                                **/
@@ -28,21 +29,22 @@ Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 /********************************************************************/
 
 GEN
-tayl(GEN x, long v, long precdl)
+tayl(GEN x, long v, long precS)
 {
-  long tetpil,i,vx = gvar9(x), av=avma;
+  long i, vx = gvar9(x);
+  gpmem_t tetpil, av=avma;
   GEN p1,y;
 
   if (v <= vx)
   {
-    long p1[] = { evaltyp(t_SER)|m_evallg(2), 0 };
-    p1[1] = evalvalp(precdl) | evalvarn(v);
+    long p1[] = { evaltyp(t_SER)|_evallg(2), 0 };
+    p1[1] = evalvalp(precS) | evalvarn(v);
     return gadd(p1,x);
   }
   p1=cgetg(v+2,t_VEC);
   for (i=0; i<v; i++) p1[i+1]=lpolx[i];
   p1[vx+1]=lpolx[v]; p1[v+1]=lpolx[vx];
-  y = tayl(changevar(x,p1), vx,precdl); tetpil=avma;
+  y = tayl(changevar(x,p1), vx,precS); tetpil=avma;
   return gerepile(av,tetpil, changevar(y,p1));
 }
 
@@ -50,17 +52,14 @@ GEN
 grando0(GEN x, long n, long do_clone)
 {
   long m, v, tx=typ(x);
-  GEN y;
 
   if (gcmp0(x)) err(talker,"zero argument in O()");
   if (tx == t_INT)
   {
     if (!gcmp1(x)) /* bug 3 + O(1). We suppose x is a truc() */
     {
-      y=cgetg(5,t_PADIC);
-      y[1] = evalvalp(n) | evalprecp(0);
-      y[2] = do_clone? lclone(x): licopy(x);
-      y[3] = un; y[4] = zero; return y;
+      if (do_clone) x = gclone(x);
+      return padiczero(x,n);
     }
     v=0; m=0; /* 1 = x^0 */
   }
@@ -68,7 +67,8 @@ grando0(GEN x, long n, long do_clone)
   {
     if (tx != t_POL && ! is_rfrac_t(tx))
       err(talker,"incorrect argument in O()");
-    v=gvar(x); m=n*gval(x,v);
+    v = gvar(x); if (v > MAXVARN) err(talker,"incorrect object in O()");
+    m = n*gval(x,v);
   }
   return zeroser(v,m);
 }
@@ -91,10 +91,12 @@ grando0(GEN x, long n, long do_clone)
 GEN
 tchebi(long n, long v) /* Assume 4*n < VERYBIGINT */
 {
-  long av,k,l;
+  long k, l;
+  gpmem_t av;
   GEN q,a,r;
 
   if (v<0) v = 0;
+  if (n < 0) err(talker,"negative degree in tchebi");
   if (n==0) return polun[v];
   if (n==1) return polx[v];
 
@@ -131,10 +133,12 @@ GEN addshiftw(GEN x, GEN y, long d);
 GEN
 legendre(long n, long v)
 {
-  long av,tetpil,m,lim;
+  long m;
+  gpmem_t av, tetpil, lim;
   GEN p0,p1,p2;
 
   if (v<0) v = 0;
+  if (n < 0) err(talker,"negative degree in legendre");
   if (n==0) return polun[v];
   if (n==1) return polx[v];
 
@@ -160,7 +164,8 @@ legendre(long n, long v)
 GEN
 cyclo(long n, long v)
 {
-  long av=avma,tetpil,d,q,m;
+  long d, q, m;
+  gpmem_t av=avma, tetpil;
   GEN yn,yd;
 
   if (n<=0) err(arither2);
@@ -249,98 +254,6 @@ roots_to_pol_r1r2(GEN a, long r1, long v)
   setlg(p1, k); return divide_conquer_prod(p1, gmul);
 }
 
-/* finds an equation for the d-th degree subfield of Q(zeta_n).
- * (Z/nZ)* must be cyclic.
- */
-GEN
-subcyclo(GEN nn, GEN dd, int v)
-{
-  long av=avma,tetpil,i,j,k,prec,q,d,p,pp,al,n,ex0,ex,aad,aa;
-  GEN a,z,pol,fa,powz,alpha;
-
-  if (typ(dd)!=t_INT || signe(dd)<=0) err(typeer,"subcyclo");
-  if (is_bigint(dd)) err(talker,"degree too large in subcyclo");
-  if (typ(nn)!=t_INT || signe(nn)<=0) err(typeer,"subcyclo");
-  if (v<0) v = 0;
-  d=itos(dd); if (d==1) return polx[v];
-  if (is_bigint(nn)) err(impl,"subcyclo for huge cyclotomic fields");
-  n = nn[2]; if ((n & 3) == 2) n >>= 1;
-  if (n == 1) err(talker,"degree does not divide phi(n) in subcyclo");
-  fa = factor(stoi(n));
-  p = itos(gmael(fa,1,1));
-  al= itos(gmael(fa,2,1));
-  if (lg((GEN)fa[1]) > 2 || (p==2 && al>2))
-    err(impl,"subcyclo in non-cyclic case");
-  if (d < n)
-  {
-    k = 1 + svaluation(d,p,&i);
-    if (k<al) { al = k; nn = gpowgs(stoi(p),al); n = nn[2]; }
-  }
-  avma=av; q = (n/p)*(p-1); /* = phi(n) */
-  if (q == d) return cyclo(n,v);
-  if (q % d) err(talker,"degree does not divide phi(n) in subcyclo");
-  q /= d;
-  if (p==2)
-  {
-    pol = powgi(polx[v],gdeux); pol[2]=un; /* replace gzero */
-    return pol; /* = x^2 + 1 */
-  }
-  a=gener(stoi(n)); aa = mael(a,2,2);
-  a=gpowgs(a,d); aad = mael(a,2,2);
-#if 1
-  prec = expi(binome(stoi(d*q-1),d)) + expi(stoi(n));
-  prec = 2 + (prec>>TWOPOTBITS_IN_LONG);
-  if (prec<DEFAULTPREC) prec = DEFAULTPREC;
-  if (DEBUGLEVEL) fprintferr("subcyclo prec = %ld\n",prec);
-  z = cgetg(3,t_COMPLEX); a=mppi(prec); setexpo(a,2); /* a = 2\pi */
-  gsincos(divrs(a,n),(GEN*)(z+2),(GEN*)(z+1),prec); /* z = e_n(1) */
-  powz = cgetg(n,t_VEC); powz[1] = (long)z;
-  k = (n+3)>>1;
-  for (i=2; i<k; i++) powz[i] = lmul(z,(GEN)powz[i-1]);
-  if ((q&1) == 0) /* totally real field, take real part */
-  {
-    for (i=1; i<k; i++) powz[i] = mael(powz,i,1);
-    for (   ; i<n; i++) powz[i] = powz[n-i];
-  }
-  else
-    for (   ; i<n; i++) powz[i] = lconj((GEN)powz[n-i]);
-
-  alpha = cgetg(d+1,t_VEC) + 1; pol=gun;
-  for (ex0=1,k=0; k<d; k++, ex0=(ex0*aa)%n)
-  {
-    GEN p1 = gzero;
-    long av1 = avma; (void)new_chunk(2*prec + 3);
-    for (ex=ex0,i=0; i<q; i++)
-    {
-      for (pp=ex,j=0; j<al; j++)
-      {
-        p1 = gadd(p1,(GEN)powz[pp]);
-        pp = mulssmod(pp,p, n);
-      }
-      ex = mulssmod(ex,aad, n);
-    }
-   /* p1 = sum z^{p^k*h}, k = 0..al-1, h runs through the subgroup of order
-    * q = phi(n)/d of (Z/nZ)^* */
-    avma = av1; alpha[k] = lneg(p1);
-  }
-  pol = roots_to_pol_intern(gun,alpha-1,v, 1);
-  if (q&1) pol=greal(pol); /* already done otherwise */
-  tetpil=avma; return gerepile(av,tetpil,ground(pol));
-#else
-{
-  /* exact computation (much slower) */
-  GEN p1 = cgetg(n+2,t_POL)+2; for (i=0; i<n; i++) p1[i]=0;
-  for (ex=1,i=0; i<q; i++, ex=(ex*aad)%n)
-    for (pp=ex,j=0; j<al; j++, pp=(pp*p)%n) p1[pp]++;
-  for (i=0; i<n; i++) p1[i] = lstoi(p1[i]);
-  p1 = normalizepol_i(p1-2,n+2); setvarn(p1,v);
-  z = cyclo(n,v); a = caract2(z,gres(p1,z),v);
-  a = gdeuc(a, modulargcd(a,derivpol(a)));
-  return gerepileupto(av, a);
-}
-#endif
-}
-
 /********************************************************************/
 /**                                                                **/
 /**                  HILBERT & PASCAL MATRICES                     **/
@@ -363,7 +276,7 @@ mathilbert(long n) /* Hilbert matrix of order n */
       coeff(p,i,j)=(long)a;
     }
   }
-  if ( n ) mael(p,1,1)=un; 
+  if ( n ) mael(p,1,1)=un;
   return p;
 }
 
@@ -371,7 +284,8 @@ mathilbert(long n) /* Hilbert matrix of order n */
 GEN
 matqpascal(long n, GEN q)
 {
-  long i,j,I, av = avma;
+  long i, j, I;
+  gpmem_t av = avma;
   GEN m, *qpow = NULL; /* gcc -Wall */
 
   if (n<0) n = -1;
@@ -411,7 +325,7 @@ matqpascal(long n, GEN q)
 GEN
 laplace(GEN x)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   long i,l,ec;
   GEN y,p1;
 
@@ -480,12 +394,9 @@ gprec(GEN x, long l)
       pr = (long) (l*pariK1+3); y=cgetr(pr); affrr(x,y); break;
 
     case t_PADIC:
-      y=cgetg(lx,tx); copyifstack(x[2], y[2]);
       if (!signe(x[4]))
-      {
-	y[1]=evalvalp(l+precp(x)) | evalprecp(0);
-	y[3]=un; y[4]=zero; return y;
-      }
+        return padiczero((GEN)x[2], l+precp(x));
+      y=cgetg(lx,tx); copyifstack(x[2], y[2]);
       y[1]=x[1]; setprecp(y,l);
       y[3]=lpuigs((GEN)x[2],l);
       y[4]=lmodii((GEN)x[4],(GEN)y[3]);
@@ -580,7 +491,8 @@ polrecip_i(GEN x)
 GEN
 binome(GEN n, long k)
 {
-  long av,i;
+  long i;
+  gpmem_t av;
   GEN y;
 
   if (k <= 1)
@@ -616,6 +528,23 @@ binome(GEN n, long k)
   return gerepileupto(av, y);
 }
 
+/* Assume n >= 1, return bin, bin[k+1] = binomial(n, k) */
+GEN
+vecbinome(long n)
+{
+  long d = (n + 1)/2, k;
+  GEN bin = cgetg(n+2, t_VEC), *C;
+  C = (GEN*)(bin + 1); /* C[k] = binomial(n, k) */
+  C[0] = gun;
+  for (k=1; k <= d; k++)
+  {
+    gpmem_t av = avma;
+    C[k] = gerepileuptoint(av, diviiexact(mulsi(n-k+1, C[k-1]), stoi(k)));
+  }
+  for (   ; k <= n; k++) C[k] = C[n - k];
+  return bin;
+}
+
 /********************************************************************/
 /**                                                                **/
 /**                  POLYNOMIAL INTERPOLATION                      **/
@@ -625,7 +554,8 @@ binome(GEN n, long k)
 GEN
 polint_i(GEN xa, GEN ya, GEN x, long n, GEN *ptdy)
 {
-  long av = avma,tetpil,i,m, ns=0, tx=typ(x);
+  long i, m, ns=0, tx=typ(x);
+  gpmem_t av = avma, tetpil;
   GEN den,ho,hp,w,y,c,d,dy;
 
   if (!xa)
@@ -675,7 +605,7 @@ GEN
 polint(GEN xa, GEN ya, GEN x, GEN *ptdy)
 {
   long tx=typ(xa), ty, lx=lg(xa);
-  
+
   if (ya) ty = typ(ya); else { ya = xa; ty = tx; xa = NULL; }
 
   if (! is_vec_t(tx) || ! is_vec_t(ty))
@@ -708,7 +638,7 @@ gtostr(GEN x)
 GEN
 gtoset(GEN x)
 {
-  ulong av;
+  gpmem_t av;
   long i,c,tx,lx;
   GEN y;
 
@@ -744,9 +674,9 @@ setisset(GEN x)
 
 /* looks if y belongs to the set x and returns the index if yes, 0 if no */
 long
-setsearch(GEN x, GEN y, long flag)
+gen_search(GEN x, GEN y, int flag, int (*cmp)(GEN,GEN))
 {
-  long av = avma,lx,j,li,ri,fl, tx = typ(x);
+  long lx,j,li,ri,fl, tx = typ(x);
 
   if (tx==t_VEC) lx = lg(x);
   else
@@ -756,22 +686,42 @@ setsearch(GEN x, GEN y, long flag)
   }
   if (lx==1) return flag? 1: 0;
 
-  if (typ(y) != t_STR) y = gtostr(y);
   li=1; ri=lx-1;
   do
   {
-    j = (ri+li)>>1; fl = gcmp((GEN)x[j],y);
-    if (!fl) { avma=av; return flag? 0: j; }
+    j = (ri+li)>>1; fl = cmp((GEN)x[j],y);
+    if (!fl) return flag? 0: j;
     if (fl<0) li=j+1; else ri=j-1;
   } while (ri>=li);
-  avma=av; if (!flag) return 0;
+  if (!flag) return 0;
   return (fl<0)? j+1: j;
 }
 
+
+long
+setsearch(GEN x, GEN y, long flag)
+{
+  gpmem_t av = avma;
+  long res;
+  if (typ(y) != t_STR) y = gtostr(y);
+  res=gen_search(x,y,flag,gcmp);
+  avma=av;
+  return res;
+}
+
+#if 0
 GEN
+gen_union(GEN x, GEN y, int (*cmp)(GEN,GEN))
+{
+  if (typ(x) != t_VEC || typ(y) != t_VEC) err(talker,"not a set in setunion");
+
+}
+#endif
+
+GEN
 setunion(GEN x, GEN y)
 {
-  long av=avma,tetpil;
+  gpmem_t av=avma, tetpil;
   GEN z;
 
   if (typ(x) != t_VEC || typ(y) != t_VEC) err(talker,"not a set in setunion");
@@ -781,7 +731,8 @@ setunion(GEN x, GEN y)
 GEN
 setintersect(GEN x, GEN y)
 {
-  long av=avma,i,lx,c;
+  long i, lx, c;
+  gpmem_t av=avma;
   GEN z;
 
   if (!setisset(x) || !setisset(y)) err(talker,"not a set in setintersect");
@@ -792,16 +743,33 @@ setintersect(GEN x, GEN y)
 }
 
 GEN
-setminus(GEN x, GEN y)
+gen_setminus(GEN set1, GEN set2, int (*cmp)(GEN,GEN))
 {
-  long av=avma,i,lx,c;
-  GEN z;
+  gpmem_t ltop=avma;
+  long find;
+  long i,j,k;
+  GEN  diff=cgetg(lg(set1),t_VEC);
+  for(i=1,j=1,k=1; i < lg(set1); i++)
+  {
+    for(find=0; j < lg(set2); j++)
+    {
+      int s=cmp((GEN)set1[i],(GEN)set2[j]);
+      if (s<0)  break ;
+      if (s>0)  continue;
+      find=1;
+    }
+    if (!find)
+      diff[k++]=set1[i];
+  }
+  setlg(diff,k);
+  return gerepilecopy(ltop,diff);
+}
 
+GEN
+setminus(GEN x, GEN y)
+{
   if (!setisset(x) || !setisset(y)) err(talker,"not a set in setminus");
-  lx=lg(x); z=cgetg(lx,t_VEC); c=1;
-  for (i=1; i<lx; i++)
-    if (setsearch(y, (GEN)x[i], 1)) z[c++] = x[i];
-  setlg(z,c); return gerepilecopy(av,z);
+  return gen_setminus(x,y,gcmp);
 }
 
 /***********************************************************************/
@@ -824,7 +792,7 @@ dirval(GEN x)
 GEN
 dirmul(GEN x, GEN y)
 {
-  ulong av = avma, lim = stack_lim(av,1);
+  gpmem_t av = avma, lim = stack_lim(av, 1);
   long lx,ly,lz,dx,dy,i,j,k;
   GEN z,p1;
 
@@ -860,7 +828,7 @@ dirmul(GEN x, GEN y)
 GEN
 dirdiv(GEN x, GEN y)
 {
-  ulong av = avma;
+  gpmem_t av = avma;
   long lx,ly,lz,dx,dy,i,j;
   GEN z,p1;
 
@@ -944,8 +912,8 @@ genrand(GEN N)
     ulong n = N[i], r;
     if (n == 0) r = 0;
     else
-    {   
-      long av = avma;
+    {
+      gpmem_t av = avma;
       if (i < nz) n++; /* allow for equality if we can go down later */
       p1 = muluu(n, gp_rand()); /* < n * 2^32, so 0 <= first word < n */
       r = (lgefint(p1)<=3)? 0: p1[2]; avma = av;
@@ -958,7 +926,7 @@ genrand(GEN N)
   i -= 2; x += i; lx -= i;
   x[1] = evalsigne(lx>2) | evallgefint(lx);
   x[0] = evaltyp(t_INT) | evallg(lx);
-  avma = (long)x; return x;
+  avma = (gpmem_t)x; return x;
 }
 
 long
@@ -982,7 +950,7 @@ gettime(void) { return timer(); }
 GEN
 numtoperm(long n, GEN x)
 {
-  ulong av;
+  gpmem_t av;
   long i,a,r;
   GEN v,w;
 
@@ -1005,7 +973,8 @@ numtoperm(long n, GEN x)
 GEN
 permtonum(GEN x)
 {
-  long av=avma, lx=lg(x)-1, n=lx, last, ind, tx = typ(x);
+  long lx=lg(x)-1, n=lx, last, ind, tx = typ(x);
+  gpmem_t av=avma;
   GEN ary,res;
 
   if (!is_vec_t(tx)) err(talker,"not a vector in permtonum");
@@ -1033,48 +1002,71 @@ permtonum(GEN x)
 /**                       MODREVERSE                               **/
 /**                                                                **/
 /********************************************************************/
+/* evaluate f(x) mod T */
+GEN
+RX_RXQ_compo(GEN f, GEN x, GEN T)
+{
+  gpmem_t av = avma;
+  long l;
+  GEN y;
 
+  if (typ(f) != t_POL) return gcopy(f);
+  l = lgef(f)-1; y = (GEN)f[l];
+  for (l--; l>=2; l--)
+    y = gres(gadd(gmul(y,x), (GEN)f[l]), T);
+  return gerepileupto(av, y);
+}
+
+/* return (1,...a^l) mod T. Not memory clean */
 GEN
+RXQ_powers(GEN a, GEN T, long l)
+{
+  long i;
+  GEN v;
+
+  if (typ(a) != t_POL) err(typeer,"RXQ_powers");
+  l += 2;
+  v = cgetg(l,t_VEC);
+  v[1] = un; if (l == 2) return v;
+
+  if (degpol(a) >= degpol(T)) a = gres(a, T);
+  v[2] = (long)a;
+  for (i=3; i<l; i++) v[i] = lres(gmul((GEN)v[i-1], a), T);
+  return v;
+}
+
+/* return y such that Mod(y, charpoly(Mod(a,T)) = Mod(a,T) */
+GEN
+modreverse_i(GEN a, GEN T)
+{
+  gpmem_t av = avma;
+  long n = degpol(T);
+  GEN y;
+
+  if (n <= 0) return gcopy(a);
+  if (n == 1)
+    return gerepileupto(av, gneg(gdiv((GEN)T[2], (GEN)T[3])));
+  if (gcmp0(a) || typ(a) != t_POL) err(talker,"reverse polmod does not exist");
+
+  y = vecpol_to_mat(RXQ_powers(a,T,n-1), n);
+  y = gauss(y, _ei(n, 2));
+  return gerepilecopy(av, vec_to_pol(y, varn(T)));
+}
+
+GEN
 polymodrecip(GEN x)
 {
-  long v,i,j,n,av,tetpil,lx;
-  GEN p1,p2,p3,p,phi,y,col;
+  long v, n;
+  GEN T, a, y;
 
-  if (typ(x)!=t_POLMOD) err(talker,"not a polymod in polymodrecip");
-  p=(GEN)x[1]; phi=(GEN)x[2];
-  v=varn(p); n=degpol(p); if (n<=0) return gcopy(x);
-  if (n==1)
-  {
-    y=cgetg(3,t_POLMOD);
-    if (typ(phi)==t_POL) phi = (GEN)phi[2];
-    p1=cgetg(4,t_POL); p1[1]=p[1]; p1[2]=lneg(phi); p1[3]=un;
-    y[1]=(long)p1;
-    if (gcmp0((GEN)p[2])) p1 = zeropol(v);
-    else
-    {
-      p1=cgetg(3,t_POL); av=avma;
-      p1[1] = evalsigne(1) | evalvarn(n) | evallgef(3);
-      p2=gdiv((GEN)p[2],(GEN)p[3]); tetpil=avma;
-      p1[2] = lpile(av,tetpil,gneg(p2));
-    }
-    y[2]=(long)p1; return y;
-  }
-  if (gcmp0(phi) || typ(phi) != t_POL)
-    err(talker,"reverse polymod does not exist");
-  av=avma; y=cgetg(n+1,t_MAT);
-  y[1]=(long)gscalcol_i(gun,n);
-  p2=phi;
-  for (j=2; j<=n; j++)
-  {
-    lx=lgef(p2); p1=cgetg(n+1,t_COL); y[j]=(long)p1;
-    for (i=1; i<=lx-2; i++) p1[i]=p2[i+1];
-    for (   ; i<=n; i++) p1[i]=zero;
-    if (j<n) p2 = gmod(gmul(p2,phi), p);
-  }
-  col=cgetg(n+1,t_COL); col[1]=zero; col[2]=un;
-  for (i=3; i<=n; i++) col[i]=zero;
-  p1=gauss(y,col); p2=gtopolyrev(p1,v); p3=caract(x,v);
-  tetpil=avma; return gerepile(av,tetpil,gmodulcp(p2,p3));
+  if (typ(x)!=t_POLMOD) err(talker,"not a polmod in modreverse");
+  T = (GEN)x[1];
+  a = (GEN)x[2];
+  n = degpol(T); if (n <= 0) return gcopy(x);
+  v = varn(T);
+  y = cgetg(3,t_POLMOD);
+  y[1] = (n==1)? lsub(polx[v], a): (long)caract2(T, a, v);
+  y[2] = (long)modreverse_i(a, T); return y;
 }
 
 /********************************************************************/
@@ -1119,7 +1111,6 @@ longcmp(GEN x, GEN y)
 
 /* Sort x = vector of elts, using cmp to compare them.
  *  flag & cmp_IND: indirect sort: return permutation that would sort x
- * For private use:
  *  flag & cmp_C  : as cmp_IND, but return permutation as vector of C-longs
  */
 GEN
@@ -1145,47 +1136,42 @@ gen_sort(GEN x, int flag, int (*cmp)(GEN,GEN))
   }
   if (!cmp) cmp = &longcmp;
   indx = (GEN) gpmalloc(lx*sizeof(long));
-  for (j=1; j<lx; j++) indx[j]=j;
+  for (j=1; j<lx; j++) indx[j] = j;
 
-  ir=lx-1; l=(ir>>1)+1;
+  ir = lx-1; l = (ir>>1)+1;
   for(;;)
   {
-    if (l>1)
-      { l--; indxt = indx[l]; }
+    if (l > 1) indxt = indx[--l];
     else
     {
-      indxt = indx[ir]; indx[ir]=indx[1]; ir--;
-      if (ir == 1)
-      {
-        indx[1] = indxt;
-        if (flag & cmp_C)
-	  for (i=1; i<lx; i++) y[i]=indx[i];
-        else if (flag & cmp_IND)
-	  for (i=1; i<lx; i++) y[i]=lstoi(indx[i]);
-        else
-	  for (i=1; i<lx; i++) y[i]=lcopy((GEN)x[indx[i]]);
-        free(indx); return y;
-      }
+      indxt = indx[ir]; indx[ir] = indx[1];
+      if (--ir == 1) break; /* DONE */
     }
-    q = (GEN)x[indxt]; i=l;
-    if (flag & cmp_REV)
-      for (j=i<<1; j<=ir; j<<=1)
-      {
-        if (j<ir && cmp((GEN)x[indx[j]],(GEN)x[indx[j+1]]) > 0) j++;
-        if (cmp(q,(GEN)x[indx[j]]) <= 0) break;
+    q = (GEN)x[indxt]; i = l;
+    for (j=i<<1; j<=ir; j<<=1)
+    {
+      if (j<ir && cmp((GEN)x[indx[j]],(GEN)x[indx[j+1]]) < 0) j++;
+      if (cmp(q,(GEN)x[indx[j]]) >= 0) break;
 
-        indx[i]=indx[j]; i=j;
-      }
-    else
-      for (j=i<<1; j<=ir; j<<=1)
-      {
-        if (j<ir && cmp((GEN)x[indx[j]],(GEN)x[indx[j+1]]) < 0) j++;
-        if (cmp(q,(GEN)x[indx[j]]) >= 0) break;
+      indx[i] = indx[j]; i = j;
+    }
+    indx[i] = indxt;
+  }
+  indx[1] = indxt;
 
-        indx[i]=indx[j]; i=j;
-      }
-    indx[i]=indxt;
+  if (flag & cmp_REV)
+  { /* reverse order */
+    GEN indy = (GEN) gpmalloc(lx*sizeof(long));
+    for (j=1; j<lx; j++) indy[j] = indx[lx-j];
+    free(indx); indx = indy;
   }
+  if (flag & cmp_C)
+    for (i=1; i<lx; i++) y[i] = indx[i];
+  else if (flag & cmp_IND)
+    for (i=1; i<lx; i++) y[i] = lstoi(indx[i]);
+  else
+    for (i=1; i<lx; i++) y[i] = lcopy((GEN)x[indx[i]]);
+  free(indx); return y;
 }
 
 #define sort_fun(flag) ((flag & cmp_LEX)? &lexcmp: &gcmp)