/********************************************************************/
/**                                                                **/
/**                  BIBLIOTHEQUE  MATHEMATIQUE                    **/
/**                     (deuxieme partie)                          **/
/**                                                                **/
/********************************************************************/
/* $Id: bibli2.c,v 1.1.1.1 1999/09/16 13:47:24 karim Exp $ */
#include "pari.h"

/********************************************************************/
/**                                                                **/
/**                 DEVELOPPEMENTS  LIMITES                        **/
/**                                                                **/
/********************************************************************/

GEN
tayl(GEN x, long v, long precdl)
{
  long tetpil,i,vx = gvar9(x), av=avma;
  GEN p1,y;

  if (v <= vx)
  {
    long p1[] = { evaltyp(t_SER)|m_evallg(2), 0 };
    p1[1] = evalvalp(precdl) | 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;
  return gerepile(av,tetpil, changevar(y,p1));
}

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;
    }
    v=0; m=0; /* 1 = x^0 */
  }
  else
  {
    if (tx != t_POL && ! is_rfrac_t(tx))
      err(talker,"incorrect argument in O()");
    v=gvar(x); m=n*gval(x,v);
  }
  return zeroser(v,m);
}

/*******************************************************************/
/**                                                               **/
/**                      SPECIAL POLYNOMIALS                      **/
/**                                                               **/
/*******************************************************************/
GEN addshiftw(GEN x, GEN y, long d);

/* Tchebichev polynomial */
/* T0=1; T1=X; T(n)=2*X*T(n-1)-T(n-2) */
GEN
tchebi(long n, long v)
{
  long av,tetpil,lim,m;
  GEN p0,p1,p2,q,z;

  if (v<0) v = 0;
  if (n==0) return polun[v];
  if (n==1) return polx[v];

  p0=gneg(polun[v]); p1=polx[v]; z = zeropol(v);
  av=avma; lim=stack_lim(av,2);
  for (m=2; m<n; m++)
  {
    q = addshiftw(gmul2n(p1,1), z, 1); tetpil=avma;
    setvarn(q,v);
    p2 = gadd(q, p0); p0=gneg(p1); p1=p2;
    if (low_stack(lim, stack_lim(av,2)))
    {
      GEN *gptr[2]; gptr[0]=&p0; gptr[1]=&p1;
      if(DEBUGMEM>1) err(warnmem,"tchebi");
      gerepilemanysp(av,tetpil,gptr,2);
    }
  }
  q = addshiftw(gmul2n(p1,1), z, 1); tetpil=avma;
  setvarn(q,v);
  return gerepile(av,tetpil,gadd(q,p0));
}

/* Legendre polynomial */
/* L0=1; L1=X; (n+1)*L(n+1)=(2*n+1)*X*L(n)-n*L(n-1) */
GEN
legendre(long n, long v)
{
  long av,tetpil,m,lim;
  GEN p0,p1,p2;

  if (v<0) v = 0;
  if (n==0) return polun[v];
  if (n==1) return polx[v];

  p0=polun[v]; av=avma; lim=stack_lim(av,2);
  p1=gmul2n(polx[v],1);
  for (m=1; m<n; m++)
  {
    p2 = addshiftw(gmulsg(4*m+2,p1), gmulsg(-4*m,p0), 1);
    setvarn(p2,v);
    p0 = p1; tetpil=avma; p1 = gdivgs(p2,m+1);
    if (low_stack(lim, stack_lim(av,2)))
    {
      GEN *gptr[2];
      if(DEBUGMEM>1) err(warnmem,"legendre");
      p0=gcopy(p0); gptr[0]=&p0; gptr[1]=&p1;
      gerepilemanysp(av,tetpil,gptr,2);
    }
  }
  tetpil=avma; return gerepile(av,tetpil,gmul2n(p1,-n));
}

/* cyclotomic polynomial */
GEN
cyclo(long n, long v)
{
  long av=avma,tetpil,d,q,m;
  GEN yn,yd;

  if (n<=0) err(arither2);
  if (v<0) v = 0;
  yn = yd = polun[0];
  for (d=1; d*d<=n; d++)
  {
    if (n%d) continue;
    q=n/d;
    m = mu(stoi(q));
    if (m)
    { /* y *= (x^d - 1) */
      if (m>0) yn = addshiftw(yn, gneg(yn), d);
      else     yd = addshiftw(yd, gneg(yd), d);
    }
    if (q==d) break;
    m = mu(stoi(d));
    if (m)
    { /* y *= (x^q - 1) */
      if (m>0) yn = addshiftw(yn, gneg(yn), q);
      else     yd = addshiftw(yd, gneg(yd), q);
    }
  }
  tetpil=avma; yn = gerepile(av,tetpil,gdeuc(yn,yd));
  setvarn(yn,v); return yn;
}

/* compute prod (L*x ± a[i]) */
GEN
roots_to_pol_intern(GEN L, GEN a, long v, int plus)
{
  long i,k,lx = lg(a), code;
  GEN p1,p2;
  if (lx == 1) return polun[v];
  p1 = cgetg(lx, t_VEC);
  code = evalsigne(1)|evalvarn(v)|evallgef(5);
  for (k=1,i=1; i<lx-1; i+=2)
  {
    p2 = cgetg(5,t_POL); p1[k++] = (long)p2;
    p2[2] = lmul((GEN)a[i],(GEN)a[i+1]);
    p2[3] = ladd((GEN)a[i],(GEN)a[i+1]);
    if (plus == 0) p2[3] = lneg((GEN)p2[3]);
    p2[4] = (long)L; p2[1] = code;
  }
  if (i < lx)
  {
    p2 = cgetg(4,t_POL); p1[k++] = (long)p2;
    p2[1] = code = evalsigne(1)|evalvarn(v)|evallgef(4);
    p2[2] = plus? a[i]: lneg((GEN)a[i]);
    p2[3] = (long)L;
  }
  setlg(p1, k); return divide_conquer_prod(p1, gmul);
}

GEN
roots_to_pol(GEN a, long v)
{
  return roots_to_pol_intern(gun,a,v,0);
}

/* prod_{i=1..r1} (x - a[i]) prod_{i=1..r2} (x - a[i])(x - conj(a[i]))*/
GEN
roots_to_pol_r1r2(GEN a, long r1, long v)
{
  long i,k,lx = lg(a), code;
  GEN p1;
  if (lx == 1) return polun[v];
  p1 = cgetg(lx, t_VEC);
  code = evalsigne(1)|evalvarn(v)|evallgef(5);
  for (k=1,i=1; i<r1; i+=2)
  {
    GEN p2 = cgetg(5,t_POL); p1[k++] = (long)p2;
    p2[2] = lmul((GEN)a[i],(GEN)a[i+1]);
    p2[3] = lneg(gadd((GEN)a[i],(GEN)a[i+1]));
    p2[4] = un; p2[1] = code;
  }
  if (i < r1+1)
    p1[k++] = ladd(polx[v], gneg((GEN)a[i]));
  for (i=r1+1; i<lx; i++)
  {
    GEN p2 = cgetg(5,t_POL); p1[k++] = (long)p2;
    p2[2] = lnorm((GEN)a[i]);
    p2[3] = lneg(gtrace((GEN)a[i]));
    p2[4] = un; p2[1] = code;
  }
  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                     **/
/**                                                                **/
/********************************************************************/
GEN addshiftpol(GEN x, GEN y, long d);

GEN
mathilbert(long n) /* Hilbert matrix of order n */
{
  long i,j;
  GEN a,p;

  if (n<0) n = 0;
  p = cgetg(n+1,t_MAT);
  for (j=1; j<=n; j++)
  {
    p[j]=lgetg(n+1,t_COL);
    for (i=1; i<=n; i++)
    {
      a=cgetg(3,t_FRAC); a[1]=un; a[2]=lstoi(i+j-1);
      coeff(p,i,j)=(long)a;
    }
  }
  return p;
}

/* q-Pascal triangle = (choose(i,j)_q) (ordinary binomial if q = NULL) */
GEN
matqpascal(long n, GEN q)
{
  long i,j,I, av = avma;
  GEN m, *qpow;

  if (n<0) n = -1;
  n++; m = cgetg(n+1,t_MAT);
  for (j=1; j<=n; j++) m[j] = lgetg(n+1,t_COL);
  if (q)
  {
    I = (n+1)/2;
    if (I > 1) { qpow = (GEN*)new_chunk(I+1); qpow[2]=q; }
    for (j=3; j<=I; j++) qpow[j] = gmul(q, qpow[j-1]);
  }
  for (i=1; i<=n; i++)
  {
    I = (i+1)/2; coeff(m,i,1)=un;
    if (q)
    {
      for (j=2; j<=I; j++)
        coeff(m,i,j) = ladd(gmul(qpow[j],gcoeff(m,i-1,j)), gcoeff(m,i-1,j-1));
    }
    else
    {
      for (j=2; j<=I; j++)
        coeff(m,i,j) = laddii(gcoeff(m,i-1,j), gcoeff(m,i-1,j-1));
    }
    for (   ; j<=i; j++) coeff(m,i,j) = coeff(m,i,i+1-j);
    for (   ; j<=n; j++) coeff(m,i,j) = zero;
  }
  return gerepileupto(av, gcopy(m));
}

/********************************************************************/
/**                                                                **/
/**                  LAPLACE TRANSFORM (OF A SERIES)               **/
/**                                                                **/
/********************************************************************/

GEN
laplace(GEN x)
{
  long i,l,ec,av,tetpil;
  GEN y,p1;

  if (typ(x)!=t_SER) err(talker,"not a series in laplace");
  if (gcmp0(x)) return gcopy(x);

  av=avma; ec=valp(x);
  if (ec<0) err(talker,"negative valuation in laplace");
  l=lg(x); y=cgetg(l,t_SER);
  p1=mpfact(ec); y[1]=x[1];
  for (i=2; i<l; i++)
  {
    y[i]=lmul(p1,(GEN)x[i]);
    ec++; p1=mulsi(ec,p1);
  }
  tetpil=avma; return gerepile(av,tetpil,gcopy(y));
}

/********************************************************************/
/**                                                                **/
/**              CONVOLUTION PRODUCT (OF TWO SERIES)               **/
/**                                                                **/
/********************************************************************/

GEN
convol(GEN x, GEN y)
{
  long l,i,j,v, vx=varn(x), lx=lg(x), ly=lg(y), ex=valp(x), ey=valp(y);
  GEN z;

  if (typ(x) != t_SER || typ(y) != t_SER)
    err(talker,"not a series in convol");
  if (gcmp0(x) || gcmp0(y))
    err(talker,"zero series in convol");
  if (varn(y) != vx)
    err(talker,"different variables in convol");
  v=ex; if (ey>v) v=ey;
  l=ex+lx; i=ey+ly; if (i<l) l=i;
  l -= v; if (l<3) err(talker,"non significant result in convol");
  for (i=v+2; i < v+l; i++)
    if (!gcmp0((GEN)x[i-ex]) && !gcmp0((GEN)y[i-ey])) { i++; break; }
  if (i == l+v) return zeroser(vx, v+l-2);

  z = cgetg(l-i+3+v,t_SER);
  z[1] = evalsigne(1) | evalvalp(i-3) | evalvarn(vx);
  for (j=i-1; j<l+v; j++) z[j-i+3]=lmul((GEN)x[j-ex],(GEN)y[j-ey]);
  return z;
}

/******************************************************************/
/**                                                              **/
/**                       PRECISION CHANGES                      **/
/**                                                              **/
/******************************************************************/

GEN
gprec(GEN x, long l)
{
  long tx=typ(x),lx=lg(x),i,pr;
  GEN y;

  if (l<=0) err(talker,"precision<=0 in gprec");
  switch(tx)
  {
    case t_REAL:
      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;
      }
      y[1]=x[1]; setprecp(y,l);
      y[3]=lpuigs((GEN)x[2],l);
      y[4]=lmodii((GEN)x[4],(GEN)y[3]);
      break;

    case t_SER:
      if (gcmp0(x)) return zeroser(varn(x), l);
      y=cgetg(l+2,t_SER); y[1]=x[1]; l++; i=l;
      if (l>=lx)
	for ( ; i>=lx; i--) y[i]=zero;
      for ( ; i>=2; i--) y[i]=lcopy((GEN)x[i]);
      break;

    case t_POL:
      lx=lgef(x); y=cgetg(lx,tx); y[1]=x[1];
      for (i=2; i<lx; i++) y[i]=lprec((GEN)x[i],l);
      break;

    case t_COMPLEX: case t_POLMOD: case t_RFRAC: case t_RFRACN:
    case t_VEC: case t_COL: case t_MAT:
      y=cgetg(lx,tx);
      for (i=1; i<lx; i++) y[i]=lprec((GEN)x[i],l);
      break;
    default: y=gcopy(x);
  }
  return y;
}

/* internal: precision given in word length (including codewords) */
GEN
gprec_w(GEN x, long pr)
{
  long tx=typ(x),lx=lg(x),i;
  GEN y;

  switch(tx)
  {
    case t_REAL:
      y=cgetr(pr); affrr(x,y); break;

    case t_POL:
      lx=lgef(x); y=cgetg(lx,tx); y[1]=x[1];
      for (i=2; i<lx; i++) y[i]=(long)gprec_w((GEN)x[i],pr);
      break;

    case t_COMPLEX: case t_POLMOD: case t_RFRAC: case t_RFRACN:
    case t_VEC: case t_COL: case t_MAT:
      y=cgetg(lx,tx);
      for (i=1; i<lx; i++) y[i]=(long)gprec_w((GEN)x[i],pr);
      break;
    default: y=gprec(x,pr);
  }
  return y;
}

/*******************************************************************/
/**                                                               **/
/**                     RECIPROCAL POLYNOMIAL                     **/
/**                                                               **/
/*******************************************************************/

GEN
polrecip(GEN x)
{
  long lx=lgef(x),i,j;
  GEN y;

  if (typ(x) != t_POL) err(typeer,"polrecip");
  y=cgetg(lx,t_POL); y[1]=x[1];
  for (i=2,j=lx-1; i<lx; i++,j--) y[i]=lcopy((GEN)x[j]);
  return normalizepol_i(y,lx);
}

/* as above. Internal (don't copy or normalize) */
GEN
polrecip_i(GEN x)
{
  long lx=lgef(x),i,j;
  GEN y;

  y=cgetg(lx,t_POL); y[1]=x[1];
  for (i=2,j=lx-1; i<lx; i++,j--) y[i]=x[j];
  return y;
}

/*******************************************************************/
/**                                                               **/
/**                      BINOMIAL COEFFICIENTS                    **/
/**                                                               **/
/*******************************************************************/

GEN
binome(GEN n, long k)
{
  long av,i;
  GEN y;

  if (k <= 1)
  {
    if (k < 0) return gzero;
    if (k == 0) return gun;
    return gcopy(n);
  }
  av = avma; y = n;
  if (typ(n) == t_INT)
  {
    if (signe(n) > 0)
    {
      GEN z = subis(n,k);
      if (cmpis(z,k) < 0) k = itos(z);
      avma = av;
      if (k <= 1)
      {
        if (k < 0) return gzero;
        if (k == 0) return gun;
        return gcopy(n);
      }
    }
    for (i=2; i<=k; i++)
      y = gdivgs(gmul(y,addis(n,i-1-k)), i);
  }
  else
  {
    for (i=2; i<=k; i++)
      y = gdivgs(gmul(y,gaddgs(n,i-1-k)), i);
  }
  return gerepileupto(av, y);
}

/********************************************************************/
/**                                                                **/
/**                  POLYNOMIAL INTERPOLATION                      **/
/**                                                                **/
/********************************************************************/

GEN
polint_i(GEN xa, GEN ya, GEN x, long n, GEN *ptdy)
{
  long av = avma,tetpil,i,m, ns=0, tx=typ(x);
  GEN den,ho,hp,w,y,c,d,dy;

  if (is_scalar_t(tx) && tx != t_INTMOD && tx != t_PADIC && tx != t_POLMOD)
  {
    GEN dif = NULL, dift;
    for (i=0; i<n; i++)
    {
      dift = gabs(gsub(x,(GEN)xa[i]), MEDDEFAULTPREC);
      if (!dif || gcmp(dift,dif)<0) { ns=i; dif=dift; }
    }
  }
  c=new_chunk(n);
  d=new_chunk(n); for (i=0; i<n; i++) c[i] = d[i] = ya[i];
  y=(GEN)d[ns--];
  for (m=1; m<n; m++)
  {
    for (i=0; i<n-m; i++)
    {
      ho = gsub((GEN)xa[i],x);
      hp = gsub((GEN)xa[i+m],x); den = gsub(ho,hp);
      if (gcmp0(den)) err(talker,"two abcissas are equal in polint");
      w=gsub((GEN)c[i+1],(GEN)d[i]); den = gdiv(w,den);
      c[i]=lmul(ho,den);
      d[i]=lmul(hp,den);
    }
    dy = (2*(ns+1) < n-m)? (GEN)c[ns+1]: (GEN)d[ns--];
    tetpil=avma; y=gadd(y,dy);
  }
  if (!ptdy) y = gerepile(av,tetpil,y);
  else
  {
    GEN *gptr[2];
    *ptdy=gcopy(dy); gptr[0]=&y; gptr[1]=ptdy;
    gerepilemanysp(av,tetpil,gptr,2);
  }
  return y;
}

GEN
polint(GEN xa, GEN ya, GEN x, GEN *ptdy)
{
  long tx=typ(xa), ty=typ(ya), lx=lg(xa);

  if (! is_vec_t(tx) || ! is_vec_t(ty))
    err(talker,"not vectors in polinterpolate");
  if (lx != lg(ya))
    err(talker,"different lengths in polinterpolate");
  if (lx <= 2)
  {
    if (lx == 1) err(talker,"no data in polinterpolate");
    ya=gcopy((GEN)ya[1]); if (ptdy) *ptdy = ya;
    return ya;
  }
  if (!x) x = polx[0];
  return polint_i(xa+1,ya+1,x,lx-1,ptdy);
}

/***********************************************************************/
/*                                                                     */
/*                          SET OPERATIONS                             */
/*                                                                     */
/***********************************************************************/

static GEN
gtostr(GEN x)
{
  char *s=GENtostr(x);
  x = strtoGENstr(s,0); free(s); return x;
}

GEN
gtoset(GEN x)
{
  long i,c,av,tetpil,tx,lx;
  GEN y;

  if (!x) return cgetg(1, t_VEC);
  tx = typ(x); lx = lg(x);
  if (!is_vec_t(tx))
  {
    if (tx != t_LIST)
      { y=cgetg(2,t_VEC); y[1]=(long)gtostr(x); return y; }
    lx = lgef(x)-1; x++;
  }
  if (lx==1) return cgetg(1,t_VEC);
  av=avma; y=cgetg(lx,t_VEC);
  for (i=1; i<lx; i++) y[i]=(long)gtostr((GEN)x[i]);
  y = sort(y);
  c=1;
  for (i=2; i<lx; i++)
    if (!gegal((GEN)y[i], (GEN)y[c])) y[++c] = y[i];
  tetpil=avma; setlg(y,c+1);
  return gerepile(av,tetpil,gcopy(y));
}

long
setisset(GEN x)
{
  long lx,i;

  if (typ(x)!=t_VEC) return 0;
  lx=lg(x)-1; if (!lx) return 1;
  for (i=1; i<lx; i++)
    if (typ(x[i]) != t_STR || gcmp((GEN)x[i+1],(GEN)x[i])<=0) return 0;
  return typ(x[i]) == t_STR;
}

/* 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)
{
  long av = avma,lx,j,li,ri,fl, tx = typ(x);

  if (tx==t_VEC) lx = lg(x);
  else
  {
    if (tx!=t_LIST) err(talker,"not a set in setsearch");
    lx=lgef(x)-1; x++;
  }
  if (lx==1) return flag? 1: 0;

  li=1; ri=lx-1;
  if (typ(y) != t_STR) y = gtostr(y);
  while (ri>=li)
  {
    j = (ri+li)>>1; fl = gcmp((GEN)x[j],y);
    if (!fl) { avma=av; return flag? 0: j; }
    if (fl<0) li=j+1; else ri=j-1;
  }
  avma=av; if (!flag) return 0;
  return (fl<0)? j+1: j;
}

GEN
setunion(GEN x, GEN y)
{
  long av=avma,tetpil;
  GEN z;

  if (typ(x) != t_VEC || typ(y) != t_VEC) err(talker,"not a set in setunion");
  z=concatsp(x,y); tetpil=avma; return gerepile(av,tetpil,gtoset(z));
}

GEN
setintersect(GEN x, GEN y)
{
  long av=avma,tetpil,i,lx,c;
  GEN z;

  if (!setisset(x) || !setisset(y)) err(talker,"not a set in setintersect");
  lx=lg(x); z=cgetg(lx,t_VEC); c=1;
  for (i=1; i<lx; i++)
    if (setsearch(y, (GEN)x[i], 0)) z[c++] = x[i];
  tetpil=avma; setlg(z,c);
  return gerepile(av,tetpil,gcopy(z));
}

GEN
setminus(GEN x, GEN y)
{
  long av=avma,tetpil,i,lx,c;
  GEN z;

  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];
  tetpil=avma; setlg(z,c);
  return gerepile(av,tetpil,gcopy(z));
}

/***********************************************************************/
/*                                                                     */
/*               OPERATIONS ON DIRICHLET SERIES                        */
/*                                                                     */
/***********************************************************************/

/* Addition, subtraction and scalar multiplication of Dirichlet series
   are done on the corresponding vectors */

static long
dirval(GEN x)
{
  long i=1,lx=lg(x);
  while (i<lx && gcmp0((GEN)x[i])) i++;
  return i;
}

GEN
dirmul(GEN x, GEN y)
{
  long lx,ly,lz,dx,dy,av,tetpil,i,j;
  GEN z,p1;

  if (typ(x)!=t_VEC || typ(y)!=t_VEC) err(talker,"not a dirseries in dirmul");
  av=avma; dx=dirval(x); dy=dirval(y); lx=lg(x); ly=lg(y);
  if (ly-dy<lx-dx) { z=y; y=x; x=z; lz=ly; ly=lx; lx=lz; lz=dy; dy=dx; dx=lz; }
  lz=min(lx*dy,ly*dx);
  z=cgetg(lz,t_VEC); for (i=1; i<lz; i++) z[i]=zero;
  for (j=dx; j<lx; j++)
  {
    p1=(GEN)x[j];
    if (!gcmp0(p1))
    {
      if (gcmp1(p1))
	for (i=j*dy; i<lz; i+=j) z[i]=ladd((GEN)z[i],(GEN)y[i/j]);
      else
      {
	if (gcmp_1(p1))
	  for (i=j*dy; i<lz; i+=j) z[i]=lsub((GEN)z[i],(GEN)y[i/j]);
	else
	  for (i=j*dy; i<lz; i+=j) z[i]=ladd((GEN)z[i],gmul(p1,(GEN)y[i/j]));
      }
    }
  }
  tetpil=avma; return gerepile(av,tetpil,gcopy(z));
}

GEN
dirdiv(GEN x, GEN y)
{
  long lx,ly,lz,dx,dy,av,tetpil,i,j;
  GEN z,p1;

  if (typ(x)!=t_VEC || typ(y)!=t_VEC) err(talker,"not a dirseries in dirmul");
  av=avma; dx=dirval(x); dy=dirval(y); lx=lg(x); ly=lg(y);
  if (dy!=1) err(talker,"not an invertible dirseries in dirdiv");
  lz=min(lx,ly*dx); p1=(GEN)y[1];
  if (!gcmp1(p1)) { y=gdiv(y,p1); x=gdiv(x,p1); }
  else x=gcopy(x);
  z=cgetg(lz,t_VEC); for (i=1; i<dx; i++) z[i]=zero;
  for (j=dx; j<lz; j++)
  {
    p1=(GEN)x[j]; z[j]=(long)p1;
    if (!gcmp0(p1))
    {
      if (gcmp1(p1))
	for (i=j+j; i<lz; i+=j) x[i]=lsub((GEN)x[i],(GEN)y[i/j]);
      else
      {
	if (gcmp_1(p1))
	  for (i=j+j; i<lz; i+=j) x[i]=ladd((GEN)x[i],(GEN)y[i/j]);
	else
	  for (i=j+j; i<lz; i+=j) x[i]=lsub((GEN)x[i],gmul(p1,(GEN)y[i/j]));
      }
    }
  }
  tetpil=avma; return gerepile(av,tetpil,gcopy(z));
}

/*************************************************************************/
/**									**/
/**			         RANDOM					**/
/**									**/
/*************************************************************************/
static long pari_randseed = 1;

/* BSD rand gives this: seed = 1103515245*seed + 12345 */
long
mymyrand()
{
#if BITS_IN_RANDOM == 64
  pari_randseed = (1000000000000654397*pari_randseed + 12347) & ~HIGHBIT;
#else
  pari_randseed = (1000276549*pari_randseed + 12347) & 0x7fffffff;
#endif
  return pari_randseed;
}

GEN muluu(ulong x, ulong y);

static ulong
gp_rand()
{
#define GLUE2(hi, lo) (((hi) << BITS_IN_HALFULONG) | (lo))
#if !defined(LONG_IS_64BIT) || BITS_IN_RANDOM == 64
  return GLUE2((mymyrand()>>12)&LOWMASK,
               (mymyrand()>>12)&LOWMASK);
#else
#define GLUE4(hi1,hi2, lo1,lo2) GLUE2(((hi1)<<16)|(hi2), ((lo1)<<16)|(lo2))
#  define LOWMASK2 0xffffUL
  return GLUE4((mymyrand()>>12)&LOWMASK2,
               (mymyrand()>>12)&LOWMASK2,
               (mymyrand()>>12)&LOWMASK2,
               (mymyrand()>>12)&LOWMASK2);
#endif
}

GEN
genrand(GEN N)
{
  long lx,i;
  GEN x;

  if (!N) return stoi(mymyrand());
  if (typ(N)!=t_INT || signe(N)<=0) err(talker,"invalid bound in random");

  lx = lgefint(N); x = new_chunk(lx);
  for (i=2; i<lx; i++)
  {
    long av = avma, n = N[i];
    ulong r = gp_rand();
    if (i < lx-1) n++; else if (!n) r = 0;
    if (n) { GEN p1 = muluu(n,r); r = (lgefint(p1)<=3)? 0: p1[2]; }
    avma = av; x[i] = r;
    if (r < (ulong)N[i]) break;
  }
  for (i++; i<lx; i++) x[i] = gp_rand();
  i=2; while (i<lx && !x[i]) i++;
  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;
}

long
setrand(long seed) { return (pari_randseed = seed); }

long
getrand() { return pari_randseed; }

long
getstack() { return top-avma; }

long
gettime() { return timer2(); }

/***********************************************************************/
/**							              **/
/**       		     PERMUTATIONS                             **/
/**								      **/
/***********************************************************************/

GEN
permute(long n, GEN x)
{
  long av=avma,i,a,r;
  GEN v,w,y;

  v=(GEN)gpmalloc((n+1)*sizeof(long)); v[1]=1;
  for (r=2; r<=n; r++)
  {
    x=dvmdis(x,r,&w); a=itos(w);
    for (i=r; i>=a+2; i--) v[i]=v[i-1];
    v[i]=r;
  }
  avma=av; y=cgetg(n+1,t_VEC);
  for (i=1; i<=n; i++) y[i]=lstoi(v[i]);
  free(v); return y;
}

GEN
permuteInv(GEN x)
{
  long av=avma,tetpil, len=lg(x)-1, ini=len, last, ind;
  GEN ary,res;

  if (typ(x)!=t_VEC && typ(x)!=t_COL) err(talker,"not a vector in permuteInv");
  res=gzero; ary=cgetg(len+1,t_VEC);
  for (ind=1; ind<=len; ind++) ary[ind]=*++x;
  ary++;
  for (last=len; last>0; last--)
  {
    len--; ind=len;
    while (ind>0 && itos((GEN)ary[ind])!=last) ind--;
    res=mulis(res,last); tetpil=avma; res=addis(res,ind);
    while (ind++<len) ary[ind-1]=ary[ind];
  }
  if (!signe(res)) { tetpil=avma; res=mpfact(ini); }
  return gerepile(av,tetpil,res);
}

/********************************************************************/
/**                                                                **/
/**                       MODREVERSE                               **/
/**                                                                **/
/********************************************************************/

GEN
polymodrecip(GEN x)
{
  long v,i,j,n,av,tetpil,lx;
  GEN p1,p2,p3,p,phi,y,col;

  if (typ(x)!=t_POLMOD) err(talker,"not a polymod in polymodrecip");
  p=(GEN)x[1]; phi=(GEN)x[2];
  v=varn(p); n=lgef(p)-3; 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));
}

/********************************************************************/
/**                                                                **/
/**                           HEAPSORT                             **/
/**                                                                **/
/********************************************************************/
static GEN vcmp_k;
static int vcmp_lk;
static int (*vcmp_cmp)(GEN,GEN);

int
pari_compare_int(int *a,int *b)
{
  return *a - *b;
}

int
pari_compare_long(long *a,long *b)
{
  return *a - *b;
}

static int
veccmp(GEN x, GEN y)
{
  int i,s;

  for (i=1; i<vcmp_lk; i++)
  {
    s = vcmp_cmp((GEN) x[vcmp_k[i]], (GEN) y[vcmp_k[i]]);
    if (s) return s;
  }
  return 0;
}

static int
longcmp(GEN x, GEN y)
{
  return ((long)x > (long)y)? 1: ((x == y)? 0: -1);
}

/* 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
gen_sort(GEN x, int flag, int (*cmp)(GEN,GEN))
{
  long i,j,indxt,ir,l,tx=typ(x),lx=lg(x);
  GEN q,y,indx;

  if (!is_matvec_t(tx) && tx != t_VECSMALL) err(typeer,"gen_sort");
  if (flag & cmp_C) tx = t_VECSMALL;
  else if (flag & cmp_IND) tx = t_VEC;
  y = cgetg(lx, tx);
  if (lx==1) return y;
  if (lx==2)
  {
    if (flag & cmp_C)
      y[1] = 1;
    else if (flag & cmp_IND)
      y[1] = un;
    else
      y[1] = lcopy((GEN)x[1]);
    return y;
  }
  if (!cmp) cmp = &longcmp;
  indx = (GEN) gpmalloc(lx*sizeof(long));
  for (j=1; j<lx; j++) indx[j]=j;

  ir=lx-1; l=(ir>>1)+1;
  for(;;)
  {
    if (l>1)
      { l--; 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;
      }
    }
    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;
    }
    indx[i]=indxt;
  }
}

#define sort_fun(flag) ((flag & cmp_LEX)? &lexcmp: &gcmp)

GEN
gen_vecsort(GEN x, GEN k, long flag)
{
  long i,j,l,t, lx = lg(x), tmp[2];

  if (lx<=2) return gen_sort(x,flag,sort_fun(flag));
  t = typ(k); vcmp_cmp = sort_fun(flag);
  if (t==t_INT)
  {
    tmp[1] = (long)k; k = tmp;
    vcmp_lk = 2;
  }
  else
  {
    if (! is_vec_t(t)) err(talker,"incorrect lextype in vecsort");
    vcmp_lk = lg(k);
  }
  l = 0;
  vcmp_k = (GEN)gpmalloc(vcmp_lk * sizeof(long));
  for (i=1; i<vcmp_lk; i++)
  {
    j = itos((GEN)k[i]);
    if (j<=0) err(talker,"negative index in vecsort");
    vcmp_k[i]=j; if (j>l) l=j;
  }
  t = typ(x);
  if (! is_matvec_t(t)) err(typeer,"vecsort");
  for (j=1; j<lx; j++)
  {
    t = typ(x[j]);
    if (! is_vec_t(t)) err(typeer,"vecsort");
    if (lg((GEN)x[j]) <= l) err(talker,"index too large in vecsort");
  }
  x = gen_sort(x, flag, veccmp);
  free(vcmp_k); return x;
}

GEN
vecsort0(GEN x, GEN k, long flag)
{
  if (flag < 0 || flag >= cmp_C) err(flagerr,"vecsort");
  return k? gen_vecsort(x,k,flag): gen_sort(x,flag, sort_fun(flag));
}

GEN
vecsort(GEN x, GEN k)
{
  return gen_vecsort(x,k, 0);
}

GEN
sindexsort(GEN x)
{
  return gen_sort(x, cmp_IND | cmp_C, gcmp);
}

GEN
sindexlexsort(GEN x)
{
  return gen_sort(x, cmp_IND | cmp_C, lexcmp);
}

GEN
indexsort(GEN x)
{
  return gen_sort(x, cmp_IND, gcmp);
}

GEN
indexlexsort(GEN x)
{
  return gen_sort(x, cmp_IND, lexcmp);
}

GEN
sort(GEN x)
{
  return gen_sort(x, 0, gcmp);
}

GEN
lexsort(GEN x)
{
  return gen_sort(x, 0, lexcmp);
}

/* index of x in table T, 0 otherwise */
long
tablesearch(GEN T, GEN x, int (*cmp)(GEN,GEN))
{
  long l=1,u=lg(T)-1,i,s;

  while (u>=l)
  {
    i = (l+u)>>1; s = cmp(x,(GEN)T[i]);
    if (!s) return i;
    if (s<0) u=i-1; else l=i+1;
  }
  return 0;
}

/* assume lg(x) = lg(y), x,y in Z^n */
int
cmp_vecint(GEN x, GEN y)
{
  long fl,i, lx = lg(x);
  for (i=1; i<lx; i++)
    if (( fl = cmpii((GEN)x[i], (GEN)y[i]) )) return fl;
  return 0;
}

/* assume x and y in primedec format. */
int
cmp_prime_over_p(GEN x, GEN y)
{
  int k = mael(x,4,2) - mael(y,4,2); /* diff. between residue degree */
  return k? ((k > 0)? 1: -1)
          : cmp_vecint((GEN)x[2], (GEN)y[2]);
}

int
cmp_prime_ideal(GEN x, GEN y)
{
  int k = cmpii((GEN)x[1], (GEN)y[1]);
  return k? k: cmp_prime_over_p(x,y);
}