/**************************************************************/
/*                                                            */
/*                        NUMBER FIELDS                       */
/*                                                            */
/**************************************************************/
/* $Id: base1.c,v 1.2 1999/09/23 17:50:55 karim Exp $ */
#include "pari.h"
#include "parinf.h"

void
checkrnf(GEN rnf)
{
  if (typ(rnf)!=t_VEC) err(idealer1);
  if (lg(rnf)!=12) err(idealer1);
}

GEN
checkbnf(GEN bnf)
{
  if (typ(bnf)!=t_VEC) err(idealer1);
  switch (lg(bnf))
  {
    case 11: return bnf;
    case 10:
      if (typ(bnf[1])==t_POL)
        err(talker,"please apply bnfinit first");
      break;
    case 7: 
      return checkbnf((GEN)bnf[1]);
      
    case 3:
      if (typ(bnf[2])==t_POLMOD)
        return checkbnf((GEN)bnf[1]);
  }
  err(idealer1);
  return NULL; /* not reached */
}

GEN
checknf(GEN nf)
{
  if (typ(nf)==t_POL) err(talker,"please apply nfinit first");
  if (typ(nf)!=t_VEC) err(idealer1);
  switch(lg(nf))
  {
    case 10: return nf;
    case 11: return checknf((GEN)nf[7]);
    case 3: if (typ(nf[2]) == t_POLMOD) return checknf((GEN)nf[1]);
  }
  err(idealer1);
  return NULL; /* not reached */
}

void
checkbnr(GEN bnr)
{
  if (typ(bnr)!=t_VEC || lg(bnr)!=7)
    err(talker,"incorrect bigray field");
  (void)checkbnf((GEN)bnr[1]);
}

void
checkbnrgen(GEN bnr)
{
  checkbnr(bnr);
  if (lg(bnr[5])<=3)
    err(talker,"please apply bnrinit(,,1) and not bnrinit(,)");
}

GEN
check_units(GEN bnf, char *f)
{
  GEN nf, x;
  bnf = checkbnf(bnf); nf = checknf(bnf);
  x = (GEN)bnf[8];
  if (lg(x) < 7 || (lg(x[5]) == 1 && lg(nf[6]) > 2))
    err(talker,"missing units in %s", f);
  return (GEN)x[5];
}

void
checkid(GEN x, long N)
{
  if (typ(x)!=t_MAT) err(idealer2);
  if (lg(x) == 1 || lg(x[1]) != N+1)
    err(talker,"incorrect matrix for ideal");
}

void
checkbid(GEN bid)
{
  if (typ(bid)!=t_VEC || lg(bid)!=6 || lg(bid[1])!=3)
    err(talker,"incorrect bigideal");
}

void
checkprimeid(GEN id)
{
  if (typ(id) != t_VEC || lg(id) != 6)
    err(talker,"incorrect prime ideal");
}

void
checkprhall(GEN prhall)
{
  if (typ(prhall) != t_VEC || lg(prhall) != 3 || typ(prhall[1]) != t_MAT)
    err(talker,"incorrect prhall format");
}

void
check_pol_int(GEN x)
{
  long k = lg(x)-1;
  for ( ; k>1; k--)
    if (typ(x[k])!=t_INT) err(talker,"polynomial not in Z[X]");
}

GEN
get_primeid(GEN x)
{
  if (typ(x) != t_VEC) return NULL;
  if (lg(x) == 3) x = (GEN)x[1];
  if (lg(x) != 6 || typ(x[3]) != t_INT) return NULL;
  return x;
}

GEN
get_bnf(GEN x, int *t)
{
  switch(typ(x))
  {
    case t_POL: *t = typ_POL;  return NULL;
    case t_QUAD: *t = typ_Q  ; return NULL;
    case t_VEC:
      switch(lg(x))
      {
        case 3:
          if (typ(x[2]) != t_POLMOD) break;
          return get_bnf((GEN)x[1],t);
        case 6 : *t = typ_QUA; return NULL;
        case 10: *t = typ_NF; return NULL;
        case 11: *t = typ_BNF; return x;
        case 7 : *t = typ_BNR;
          x = (GEN)x[1]; if (typ(x)!=t_VEC || lg(x)!=11) break;
          return x;
      }
    case t_MAT:
      if (lg(x)==2)
        switch(lg(x[1]))
        {
          case 8: case 11:
            *t = typ_CLA; return NULL;
        }
  }
  *t = typ_NULL; return NULL;
}

GEN
get_nf(GEN x, int *t)
{
  switch(typ(x))
  {
    case t_POL : *t = typ_POL; return NULL;
    case t_QUAD: *t = typ_Q  ; return NULL;
    case t_VEC:
      switch(lg(x))
      {
        case 3:
          if (typ(x[2]) != t_POLMOD) break;
          return get_nf((GEN)x[1],t);
        case 10: *t = typ_NF; return x;
        case 11: *t = typ_BNF;
          x = (GEN)x[7]; if (typ(x)!=t_VEC || lg(x)!=10) break;
          return x;
        case 7 : *t = typ_BNR;
          x = (GEN)x[1]; if (typ(x)!=t_VEC || lg(x)!=11) break;
          x = (GEN)x[7]; if (typ(x)!=t_VEC || lg(x)!=10) break;
          return x;

        case 14: case 20:
          *t = typ_ELL; return NULL;
      }break;
    case t_MAT:
      if (lg(x)==2)
        switch(lg(x[1]))
        {
          case 8: case 11:
            *t = typ_CLA; return NULL;
        }
  }
  *t = typ_NULL; return NULL;
}

/*************************************************************************/
/**									**/
/**			       GALOIS GROUP   				**/
/**									**/
/*************************************************************************/

/* exchange elements i and j in vector x */
static GEN
transroot(GEN x, int i, int j)
{
  long k;
  x = dummycopy(x);
  k=x[i]; x[i]=x[j]; x[j]=k; return x;
}

#define randshift (BITS_IN_RANDOM - 3)

GEN
tschirnhaus(GEN x)
{
  long a, v = varn(x), av = avma;
  GEN u, p1 = cgetg(5,t_POL);

  if (typ(x)!=t_POL) err(notpoler,"tschirnhaus");
  if (lgef(x) < 4) err(constpoler,"tschirnhaus");
  if (v) { u=dummycopy(x); setvarn(u,0); x=u; }
  p1[1] = evalsigne(1)|evalvarn(0)|evallgef(5);
  do
  {
    a = mymyrand() >> randshift; if (a==0) a =1; p1[4]=lstoi(a);
    a = mymyrand() >> (randshift-1); if (a>=4) a-=8; p1[3]=lstoi(a);
    a = mymyrand() >> (randshift-1); if (a>=4) a-=8; p1[2]=lstoi(a);
    u = caract2(x,p1,v); a=avma;
  }
  while (lgef(srgcd(u,derivpol(u))) > 3); /* while u not separable */
  if (DEBUGLEVEL>1)
    fprintferr("Tschirnhaus transform. New pol: %Z",u);
  avma=a; return gerepileupto(av,u);
}
#undef randshift

int
gpolcomp(GEN p1, GEN p2)
{
  int s,j = lgef(p1)-2;

  if (lgef(p2)-2 != j)
    err(bugparier,"gpolcomp (different degrees)");
  for (; j>=2; j--)
  {
    s = absi_cmp((GEN)p1[j], (GEN)p2[j]);
    if (s) return s;
  }
  return 0;
}

/* pol is assumed to be non-zero, primitive and integral */
static GEN
primitive_pol_to_monic(GEN pol, GEN *ptlead)
{
  long n = lgef(pol)-1;
  GEN p1,lead,res;

  if (gcmp1((GEN)pol[n])) { if (ptlead) *ptlead = NULL; return pol; }
  res = cgetg(n+1,t_POL); res[1]=pol[1];
  lead = (GEN) pol[n]; p1 = gun;
  res[n] = un; n--;
  if (n>1) res[n] = pol[n];
  for (n--; n>1; n--)
  {
    p1 = mulii(lead,p1);
    res[n] = lmulii(p1,(GEN)pol[n]);
  }
  if (ptlead) *ptlead = lead; return res;
}

/* compute x1*x2^2 + x2*x3^2 + x3*x4^2 + x4*x1^2 */
static GEN
get_F4(GEN x)
{
  GEN p1=gzero;
  long i;

  for (i=1; i<=4; i++)
    p1 = gadd(p1, gmul((GEN)x[i], gsqr((GEN)x[(i&3)+1])));
  return p1;
}

GEN galoisbig(GEN x, long prec);
GEN roots_to_pol(GEN a, long v);

GEN
galois(GEN x, long prec)
{
  long av=avma,av1,i,j,k,n,f,l,l2,e,e1,pr,ind;
  GEN x1,p1,p2,p3,p4,p5,p6,w,y,z,ee;
  static int ind5[20]={2,5,3,4, 1,3,4,5, 1,5,2,4, 1,2,3,5, 1,4,2,3};
  static int ind6[60]={3,5,4,6, 2,6,4,5, 2,3,5,6, 2,4,3,6, 2,5,3,4,
                       1,4,5,6, 1,5,3,6, 1,6,3,4, 1,3,4,5, 1,6,2,5,
                       1,2,4,6, 1,5,2,4, 1,3,2,6, 1,2,3,5, 1,4,2,3};
  if (typ(x)!=t_POL) err(notpoler,"galois");
  n=lgef(x)-3; if (n<=0) err(constpoler,"galois");
  if (n>11) err(impl,"galois of degree higher than 11");
  x = gdiv(x,content(x));
  for (i=2; i<=n+2; i++)
    if (typ(x[i])!=t_INT) err(polrationer,"galois");
  if (gisirreducible(x) != gun)
    err(impl,"galois of reducible polynomial");

  if (n<4)
  {
    if (n<3)
    {
      avma=av; y=cgetg(4,t_VEC);
      y[1] = (n==1)? un: deux;
      y[2]=lnegi(gun);
    }
    else /* n=3 */
    {
      f=carreparfait(discsr(x));
      avma=av; y=cgetg(4,t_VEC);
      if (f) { y[1]=lstoi(3); y[2]=un; }
      else   { y[1]=lstoi(6); y[2]=lnegi(gun); }
    }
    y[3]=un; return y;
  }
  x1 = x = primitive_pol_to_monic(x,NULL); av1=avma;
  if (n>7) return galoisbig(x,prec);
  for(;;)
  {
    switch(n)
    {
      case 4: z = cgetg(7,t_VEC);
        for(;;)
	{
	  p1=roots(x,prec);
          z[1] = (long)get_F4(p1);
          z[2] = (long)get_F4(transroot(p1,1,2));
          z[3] = (long)get_F4(transroot(p1,1,3));
          z[4] = (long)get_F4(transroot(p1,1,4));
          z[5] = (long)get_F4(transroot(p1,2,3));
          z[6] = (long)get_F4(transroot(p1,3,4));
          p4 = roots_to_pol(z,0);
          p5=grndtoi(greal(p4),&e);
          e1=gexpo(gimag(p4)); if (e1>e) e=e1;
          if (e <= -10) break;
	  prec = (prec<<1)-2;
	}
	p6 = modulargcd(p5, derivpol(p5));
	if (typ(p6)==t_POL && lgef(p6)>3) goto tchi;
	p2 = (GEN)factor(p5)[1];
	switch(lg(p2)-1)
	{
	  case 1: f=carreparfait(discsr(x)); avma=av; y=cgetg(4,t_VEC);
	    y[3]=un;
	    if (f) { y[2]=un; y[1]=lstoi(12); return y; }
	    y[2]=lnegi(gun); y[1]=lstoi(24); return y;

	  case 2: avma=av; y=cgetg(4,t_VEC);
	    y[3]=un; y[2]=lnegi(gun); y[1]=lstoi(8); return y;
	
	  case 3: avma=av; y=cgetg(4,t_VEC);
	    y[1]=lstoi(4); y[3]=un;
	    y[2] = (lgef(p2[1])==5)? un: lnegi(gun);
	    return y;

	  default: err(bugparier,"galois (bug1)");
	}

      case 5: z = cgetg(7,t_VEC);
        ee = new_chunk(7); w = new_chunk(7);
        for(;;)
	{
          for(;;)
	  {
	    p1=roots(x,prec);
	    for (l=1; l<=6; l++)
	    {
	      p2=(l==1)?p1: ((l<6)?transroot(p1,1,l): transroot(p1,2,5));
	      p3=gzero;
              for (k=0,i=1; i<=5; i++,k+=4)
	      {
		p5 = gadd(gmul((GEN)p2[ind5[k]],(GEN)p2[ind5[k+1]]),
		          gmul((GEN)p2[ind5[k+2]],(GEN)p2[ind5[k+3]]));
		p3 = gadd(p3, gmul(gsqr((GEN)p2[i]),p5));
	      }
	      w[l]=lrndtoi(greal(p3),&e);
              e1 = gexpo(gimag(p3)); if (e1>e) e=e1;
              ee[l]=e; z[l] = (long)p3;
	    }
            p4 = roots_to_pol(z,0);
	    p5=grndtoi(greal(p4),&e);
            e1 = gexpo(gimag(p4)); if (e1>e) e=e1;
            if (e <= -10) break;
	    prec = (prec<<1)-2;
	  }
	  p6 = modulargcd(p5,derivpol(p5));
	  if (typ(p6)==t_POL && lgef(p6)>3) goto tchi;
	  p3=(GEN)factor(p5)[1];
	  f=carreparfait(discsr(x));
	  if (lg(p3)-1==1)
	  {
	    avma=av; y=cgetg(4,t_VEC); y[3]=un;
	    if (f) { y[2]=un; y[1]=lstoi(60); return y; }
	    else { y[2]=lneg(gun); y[1]=lstoi(120); return y; }
	  }
	  if (!f)
	  {
	    avma=av; y=cgetg(4,t_VEC);
	    y[3]=un; y[2]=lneg(gun); y[1]=lstoi(20); return y;
	  }
          pr = - (bit_accuracy(prec) >> 1);
          for (l=1; l<=6; l++)
	    if (ee[l] <= pr && gcmp0(poleval(p5,(GEN)w[l]))) break;
	  if (l>6) err(bugparier,"galois (bug4)");
	  p2=(l==6)? transroot(p1,2,5):transroot(p1,1,l);
	  p3=gzero;
	  for (i=1; i<=5; i++)
	  {
	    j=(i%5)+1;
	    p3=gadd(p3,gmul(gmul((GEN)p2[i],(GEN)p2[j]),
			    gsub((GEN)p2[j],(GEN)p2[i])));
	  }
	  p5=gsqr(p3); p4=grndtoi(greal(p5),&e);
          e1 = gexpo(gimag(p5)); if (e1>e) e=e1;
	  if (e <= -10)
	  {
	    if (gcmp0(p4)) goto tchi;
	    f=carreparfait(p4); avma=av; y=cgetg(4,t_VEC);
	    y[3]=y[2]=un; y[1]=lstoi(f?5:10);
	    return y;
	  }
	  prec=(prec<<1)-2;
	}

      case 6: z = cgetg(7, t_VEC);
        for(;;)
	{
          for(;;)
	  {
	    p1=roots(x,prec);
	    for (l=1; l<=6; l++)
	    {
	      p2=(l==1)?p1:transroot(p1,1,l);
	      p3=gzero; k=0;
              for (i=1; i<=5; i++) for (j=i+1; j<=6; j++)
	      {
		p5=gadd(gmul((GEN)p2[ind6[k]],(GEN)p2[ind6[k+1]]),
		        gmul((GEN)p2[ind6[k+2]],(GEN)p2[ind6[k+3]]));
		p3=gadd(p3,gmul(gsqr(gmul((GEN)p2[i],(GEN)p2[j])),p5)); k+=4;
	      }
	      z[l] = (long)p3;
	    }
            p4=roots_to_pol(z,0);
	    p5=grndtoi(greal(p4),&e);
            e1 = gexpo(gimag(p4)); if (e1>e) e=e1;
            if (e <= -10) break;
	    prec=(prec<<1)-2;
	  }
	  p6 = modulargcd(p5,derivpol(p5));
	  if (typ(p6)==t_POL && lgef(p6)>3) goto tchi;
	  p2=(GEN)factor(p5)[1];
	  switch(lg(p2)-1)
	  {
	    case 1:
              z = cgetg(11,t_VEC); ind=0;
	      p3=gadd(gmul(gmul((GEN)p1[1],(GEN)p1[2]),(GEN)p1[3]),
	              gmul(gmul((GEN)p1[4],(GEN)p1[5]),(GEN)p1[6]));
              z[++ind] = (long)p3;
	      for (i=1; i<=3; i++)
		for (j=4; j<=6; j++)
		{
		  p2=transroot(p1,i,j);
		  p3=gadd(gmul(gmul((GEN)p2[1],(GEN)p2[2]),(GEN)p2[3]),
		          gmul(gmul((GEN)p2[4],(GEN)p2[5]),(GEN)p2[6]));
		  z[++ind] = (long)p3;
		}
              p4 = roots_to_pol(z,0);
	      p5=grndtoi(greal(p4),&e);
              e1 = gexpo(gimag(p4)); if (e1>e) e=e1;
	      if (e <= -10)
	      {
		p6 = modulargcd(p5,derivpol(p5));
		if (typ(p6)==t_POL && lgef(p6)>3) goto tchi;
		p2=(GEN)factor(p5)[1];
		f=carreparfait(discsr(x));
		avma=av; y=cgetg(4,t_VEC); y[3]=un;
		if (lg(p2)-1==1)
		{
		  if (f) { y[2]=un; y[1]=lstoi(360); }
		  else { y[2]=lnegi(gun); y[1]=lstoi(720); }
		}
		else
		{
		  if (f) { y[2]=un; y[1]=lstoi(36); }
		  else { y[2]=lnegi(gun); y[1]=lstoi(72); }
		}
                return y;
	      }
	      prec=(prec<<1)-2; break;
		
	    case 2: l2=lgef(p2[1])-3; if (l2>3) l2=6-l2;
	      switch(l2)
	      {
		case 1: f=carreparfait(discsr(x));
		  avma=av; y=cgetg(4,t_VEC); y[3]=un;
		  if (f) { y[2]=un; y[1]=lstoi(60); }
		  else { y[2]=lneg(gun); y[1]=lstoi(120); }
		  return y;
		case 2: f=carreparfait(discsr(x));
		  if (f)
		  {
		    avma=av; y=cgetg(4,t_VEC);
		    y[3]=y[2]=un; y[1]=lstoi(24);
		  }
		  else
		  {
		    p3=(lgef(p2[1])==5) ? (GEN)p2[2]:(GEN)p2[1];
		    f=carreparfait(discsr(p3));
		    avma=av; y=cgetg(4,t_VEC); y[2]=lneg(gun);
		    if (f) { y[1]=lstoi(24); y[3]=deux; }
		    else { y[1]=lstoi(48); y[3]=un; }
		  }
		  return y;
		case 3: f=carreparfait(discsr((GEN)p2[1]))
		       || carreparfait(discsr((GEN)p2[2]));
		  avma=av; y=cgetg(4,t_VEC);
		  y[3]=un; y[2]=lneg(gun); y[1]=lstoi(f? 18: 36);
		  return y;
	      }
	    case 3:
	      for (l2=1; l2<=3; l2++)
		if (lgef(p2[l2])>=6) p3=(GEN)p2[l2];
	      if (lgef(p3)==6)
	      {
		f=carreparfait(discsr(p3)); avma=av; y=cgetg(4,t_VEC);
                y[2]=lneg(gun); y[1]=lstoi(f? 6: 12);
	      }
	      else
	      {
		f=carreparfait(discsr(x)); avma=av; y=cgetg(4,t_VEC);
		if (f) { y[2]=un; y[1]=lstoi(12); }
		else { y[2]=lneg(gun); y[1]=lstoi(24); }
	      }
              y[3]=un; return y;
	    case 4: avma=av; y=cgetg(4,t_VEC);
	      y[1]=lstoi(6); y[2]=lneg(gun); y[3]=deux; return y;
            default: err(bugparier,"galois (bug3)");
	  }
	}
	
      case 7: z = cgetg(36,t_VEC);
        for(;;)
	{
	  ind = 0; p1=roots(x,prec); p4=gun;
	  for (i=1; i<=5; i++)
	    for (j=i+1; j<=6; j++)
            {
              p6 = gadd((GEN)p1[i],(GEN)p1[j]);
	      for (k=j+1; k<=7; k++)
                z[++ind] = (long) gadd(p6,(GEN)p1[k]);
            }
          p4 = roots_to_pol(z,0);
          p5 = grndtoi(greal(p4),&e);
          e1 = gexpo(gimag(p4)); if (e1>e) e=e1;
	  if (e <= -10) break;
          prec = (prec<<1)-2;
	}
	p6 = modulargcd(p5,derivpol(p5));
	if (typ(p6)==t_POL && lgef(p6)>3) goto tchi;
	p2=(GEN)factor(p5)[1];
	switch(lg(p2)-1)
	{
	  case 1: f=carreparfait(discsr(x)); avma=av; y=cgetg(4,t_VEC); y[3]=un;
	    if (f) { y[2]=un; y[1]=lstoi(2520); }
	    else { y[2]=lneg(gun); y[1]=lstoi(5040); }
	    return y;
	  case 2: f=lgef(p2[1])-3; avma=av; y=cgetg(4,t_VEC); y[3]=un;
	    if (f==7 || f==28) { y[2]=un; y[1]=lstoi(168); }
	    else { y[2]=lneg(gun); y[1]=lstoi(42); }
	    return y;
	  case 3: avma=av; y=cgetg(4,t_VEC);
	    y[3]=y[2]=un; y[1]=lstoi(21); return y;
	  case 4: avma=av; y=cgetg(4,t_VEC);
	    y[3]=un; y[2]=lneg(gun); y[1]=lstoi(14); return y;
	  case 5: avma=av; y=cgetg(4,t_VEC);
	    y[3]=y[2]=un; y[1]=lstoi(7); return y;
          default: err(talker,"galois (bug2)");
	}
    }
    tchi: avma=av1; x=tschirnhaus(x1);
  }
}

GEN
galoisapply(GEN nf, GEN aut, GEN x)
{
  long av=avma,tetpil,lx,j,N;
  GEN p,p1,y,pol;

  nf=checknf(nf); pol=(GEN)nf[1];
  if (typ(aut)==t_POL) aut = gmodulcp(aut,pol);
  else
  {
    if (typ(aut)!=t_POLMOD || !gegal((GEN)aut[1],pol))
      err(talker,"incorrect galois automorphism in galoisapply");
  }
  switch(typ(x))
  {
    case t_INT: case t_INTMOD: case t_FRAC: case t_FRACN: case t_PADIC:
      avma=av; return gcopy(x);

    case t_POLMOD: x = (GEN) x[2]; /* fall through */
    case t_POL:
      p1 = gsubst(x,varn(pol),aut);
      if (typ(p1)!=t_POLMOD || !gegal((GEN)p1[1],pol)) p1 = gmodulcp(p1,pol);
      return gerepileupto(av,p1);

    case t_VEC:
      if (lg(x)==3)
      {
	y=cgetg(3,t_VEC);
	y[1]=(long)galoisapply(nf,aut,(GEN)x[1]);
        y[2]=lcopy((GEN)x[2]); return gerepileupto(av, y);
      }
      if (lg(x)!=6) err(typeer,"galoisapply");
      y=cgetg(6,t_VEC); y[1]=x[1]; y[3]=x[3]; y[4]=x[4];
      p = (GEN)x[1];
      p1=centermod(galoisapply(nf,aut,(GEN)x[2]), p);
      if (is_pm1(x[3]))
	if (ggval(subres(gmul((GEN)nf[7],p1),pol), p) > itos((GEN)x[4]))
	  p1[1] =  (signe(p1[1]) > 0)? lsub((GEN)p1[1], p)
	                             : ladd((GEN)p1[1], p);
      y[2]=(long)p1;
      y[5]=(long)centermod(galoisapply(nf,aut,(GEN)x[5]), p);
      tetpil=avma; return gerepile(av,tetpil,gcopy(y));

    case t_COL:
      N=lgef(pol)-3;
      if (lg(x)!=N+1) err(typeer,"galoisapply");
      p1=galoisapply(nf,aut,gmul((GEN)nf[7],x)); tetpil=avma;
      return gerepile(av,tetpil,algtobasis(nf,p1));

    case t_MAT:
      lx=lg(x); if (lx==1) return cgetg(1,t_MAT);
      N=lgef(pol)-3;
      if (lg(x[1])!=N+1) err(typeer,"galoisapply");
      p1=cgetg(lx,t_MAT);
      for (j=1; j<lx; j++) p1[j]=(long)galoisapply(nf,aut,(GEN)x[j]);
      if (lx==N+1) p1 = idealhermite(nf,p1);
      return gerepileupto(av,p1);
  }
  err(typeer,"galoisapply");
  return NULL; /* not reached */
}

GEN pol_to_monic(GEN pol, GEN *lead);
int cmp_pol(GEN x, GEN y);

static GEN
get_nfpol(GEN x, GEN *nf)
{
  if (typ(x) == t_POL) { *nf = NULL; return x; }
  *nf = checknf(x); return (GEN)(*nf)[1];
}

/* if fliso test for isomorphism, for inclusion otherwise. */
static GEN
nfiso0(GEN a, GEN b, long fliso)
{
  long av=avma,tetpil,n,m,i,vb,lx;
  GEN nfa,nfb,p1,y,la,lb;

  a = get_nfpol(a, &nfa);
  b = get_nfpol(b, &nfb);
  if (fliso && nfa && !nfb) { p1=a; a=b; b=p1; p1=nfa; nfa=nfb; nfb=p1; }
  m=lgef(a)-3;
  n=lgef(b)-3; if (m<=0 || n<=0) err(constpoler,"nfiso or nfincl");
  if (fliso)
    { if (n!=m) return gzero; }
  else
    { if (n%m) return gzero; }

  if (nfb) lb = NULL; else b = pol_to_monic(b,&lb);
  if (nfa) la = NULL; else a = pol_to_monic(a,&la);
  if (nfa && nfb)
  {
    if (fliso)
    {
      if (!gegal((GEN)nfa[2],(GEN)nfb[2])
       || !gegal((GEN)nfa[3],(GEN)nfb[3])) return gzero;
    }
    else
      if (!divise((GEN)nfb[3], gpowgs((GEN)nfa[3],n/m))) return gzero;
  }
  else
  {
    GEN da = nfa? (GEN)nfa[3]: discsr(a);
    GEN db = nfb? (GEN)nfb[3]: discsr(b);
    if (fliso)
    {
      p1=gdiv(da,db);
      if (typ(p1)==t_FRAC) p1=mulii((GEN)p1[1],(GEN)p1[2]);
      if (!carreparfait(p1)) { avma=av; return gzero; }
    }
    else
    {
      long q=n/m;
      GEN fa=factor(da), ex=(GEN)fa[2];
      fa=(GEN)fa[1]; lx=lg(fa);
      for (i=1; i<lx; i++)
        if (mod2((GEN)ex[i]) && !divise(db,gpowgs((GEN)fa[i],q)))
          { avma=av; return gzero; }
    }
  }
  a = dummycopy(a); setvarn(a,0);
  b = dummycopy(b); vb=varn(b);
  if (nfb)
  {
    if (vb == 0) nfb = gsubst(nfb, 0, polx[MAXVARN]);
    y = lift_intern(nfroots(nfb,a));
  }
  else
  {
    if (vb == 0) setvarn(b, fetch_var());
    y = (GEN)polfnf(a,b)[1]; lx = lg(y);
    for (i=1; i<lx; i++)
    {
      if (lgef(y[i]) != 4) { setlg(y,i); break; }
      y[i] = (long)gneg_i(lift_intern(gmael(y,i,2)));
    }
    y = gen_sort(y, 0, cmp_pol);
    settyp(y, t_VEC);
    if (vb == 0) delete_var();
  }
  lx = lg(y); if (lx==1) { avma=av; return gzero; }
  for (i=1; i<lx; i++)
  {
    p1 = (GEN)y[i]; setvarn(p1, vb);
    if (lb) p1 = poleval(p1, gmul(polx[vb],lb));
    y[i] = la? ldiv(p1,la): (long)p1;
  }
  tetpil=avma; return gerepile(av,tetpil,gcopy(y));
}

GEN
nfisisom(GEN a, GEN b)
{
  return nfiso0(a,b,1);
}

GEN
nfisincl(GEN a, GEN b)
{
  return nfiso0(a,b,0);
}

/*************************************************************************/
/**									**/
/**			       INITALG					**/
/**									**/
/*************************************************************************/
GEN LLL_nfbasis(GEN *x, GEN polr, GEN base, long prec);
GEN eleval(GEN f,GEN h,GEN a);
int canon_pol(GEN z);
GEN mat_to_vecpol(GEN x, long v);
GEN vecpol_to_mat(GEN v, long n);

/* For internal use. compute trace(x mod pol), sym=polsym(pol,deg(pol)-1) */
GEN
quicktrace(GEN x, GEN sym)
{
  GEN p1 = gzero;
  long i;

  if (signe(x))
  {
    sym--;
    for (i=lgef(x)-1; i>1; i--)
      p1 = gadd(p1, gmul((GEN)x[i],(GEN)sym[i]));
  }
  return p1;
}

/* Seek a new, simpler, polynomial pol defining the same number field as
 * x (assumed to be monic at this point).
 * *ptx   receives pol
 * *ptdx  receives disc(pol)
 * *ptrev expresses the new root in terms of the old one.
 * base updated in place.
 * prec = 0 <==> field is totally real
 */
static void
nfinit_reduce(long flag, GEN *px, GEN *pdx, GEN *prev, GEN *pbase, long prec)
{
  GEN chbas,a,phimax,dxn,s,sn,p1,p2,p3,polmax,rev,polr;
  GEN x = *px, dx = *pdx, base = *pbase;
  long i,j,nmax,numb,flc,v=varn(x), n=lg(base)-1;

  if (n == 1)
  {
    *px = gsub(polx[v],gun); *pdx = gun;
    *prev = polymodrecip(gmodulcp(gun, x)); return;
  }
  polr = prec? roots(x, prec): NULL;
  if (polr)
    for (s=gzero,i=1; i<=n; i++) s = gadd(s,gnorm((GEN)polr[i]));
  else
    s = subii(sqri((GEN)x[n+1]), shifti((GEN)x[n],1));
    
  chbas = LLL_nfbasis(&x,polr,base,prec);
  if (DEBUGLEVEL) msgtimer("LLL basis");

  phimax=polx[v]; polmax=dummycopy(x);
  nmax=(flag & nf_PARTIAL)?min(n,3):n;

  for (numb=0,i=2; i<=nmax || (!numb && i<=n); i++)
  { /* cf pols_for_polred */
    if (DEBUGLEVEL>2) { fprintferr("i = %ld\n",i); flusherr(); }
    p1 = a = gmul(base,(GEN)chbas[i]); p3=content(p1);
    if (gcmp1(p3)) p3 = NULL; else p1 = gdiv(p1,p3);
    p1 = caract2(x,p1,v);
    if (p3)
      for (p2=gun, j=lgef(p1)-2; j>=2; j--)
      {
        p2 = gmul(p2,p3); p1[j] = lmul((GEN)p1[j], p2);
      }
    p2 = modulargcd(derivpol(p1),p1);
    if (lgef(p2) > 3) continue;

    if (DEBUGLEVEL>3) outerr(p1);
    dxn=discsr(p1); flc=absi_cmp(dxn,dx); numb++;
    if (flc>0) continue;

    if (polr)
      for (sn=gzero,j=1; j<=n; j++)
        sn = gadd(sn,gnorm(poleval(a,(GEN)polr[j])));
    else
      sn = subii(sqri((GEN)p1[n+1]), shifti((GEN)p1[n],1));
    if (flc>=0)
    {
      flc=gcmp(sn,s);
      if (flc>0 || (!flc && gpolcomp(p1,polmax) >= 0)) continue;
    }
    dx=dxn; s=sn; polmax=p1; phimax=a;
  }
  if (!numb) 
  {
    if (gisirreducible(x)!=gun) err(redpoler,"nfinit_reduce");
    err(talker,"you have found a counter-example to a conjecture, "
               "please send us\nthe polynomial as soon as possible");
  }
  if (phimax == polx[v]) /* no improvement */
    rev = gmodulcp(phimax,x);
  else
  {
    if (canon_pol(polmax) < 0) phimax = gneg_i(phimax);
    if (DEBUGLEVEL>1) fprintferr("polmax = %Z\n",polmax);
    p2 = gmodulcp(phimax,x); rev = polymodrecip(p2);
    a = base; base = cgetg(n+1,t_VEC);
    p1 = (GEN)rev[2];
    for (i=1; i<=n; i++)
      base[i] = (long)eleval(polmax, (GEN)a[i], p1);
    p1 = vecpol_to_mat(base,n); p2=denom(p1); p1=gmul(p2,p1);
    p1 = gdiv(hnfmod(p1,detint(p1)), p2);
    base = mat_to_vecpol(p1,v);
  }
  *px=polmax; *pdx=dx; *prev=rev; *pbase = base;
}

/* pol belonging to Z[x], return a monic polynomial generating the same field
 * as pol (x-> ax+b)) set lead = NULL if pol was monic (after dividing
 * by the content), and to to leading coeff otherwise.
 * No garbage collecting done.
 */
GEN
pol_to_monic(GEN pol, GEN *lead)
{
  long n = lgef(pol)-1;
  GEN p1;

  if (n==1 || gcmp1((GEN)pol[n])) { *lead = NULL; return pol; }

  p1=content(pol); if (!gcmp1(p1)) pol = gdiv(pol,p1);
  return primitive_pol_to_monic(pol,lead);
}

GEN
make_TI(GEN nf, GEN TI, GEN con)
{
  GEN z, p1 = hnfmod(TI,detint(TI)), d = dethnf(p1);
  long n = lg(TI)-1;

  d = mulii(d, gpowgs(con, n));
  p1 = ideal_two_elt(nf, p1); p1 = gmul(p1,con);
  z = cgetg(4,t_VEC);
  z[1]=p1[1]; z[2]=p1[2]; z[3]=(long)d; return z;
}

/* basis = integer basis. roo = real part of the roots */
GEN
make_M(long n,long ru,GEN basis,GEN roo)
{
  GEN p1,res = cgetg(n+1,t_MAT);
  long i,j;

  for (j=1; j<=n; j++)
  {
    p1=cgetg(ru+1,t_COL); res[j]=(long)p1;
    for (i=1; i<=ru; i++)
      p1[i]=(long)poleval((GEN)basis[j],(GEN)roo[i]);
  }
  if (DEBUGLEVEL>4) msgtimer("matrix M");
  return res;
}

GEN
make_MC(long n,long r1,long ru,GEN M)
{
  GEN p1,p2,res=cgetg(ru+1,t_MAT);
  long i,j,av,tetpil;

  for (j=1; j<=ru; j++)
  {
    p1=cgetg(n+1,t_COL); res[j]=(long)p1;
    for (i=1; i<=n; i++)
    {
      av=avma; p2=gconj(gcoeff(M,j,i)); tetpil=avma;
      p1[i] = (j<=r1)? (long)p2: lpile(av,tetpil,gmul2n(p2,1));
    }
  }
  if (DEBUGLEVEL>4) msgtimer("matrix MC");
  return res;
}

GEN
get_roots(GEN x,long r1,long ru,long prec)
{
  GEN roo = (typ(x)==t_VEC)? dummycopy(x): roots(x,prec);
  long i;

  for (i=1; i<=r1; i++) roo[i]=lreal((GEN)roo[i]);
  for (   ; i<=ru; i++) roo[i]=roo[(i<<1)-r1];
  roo[0]=evaltyp(t_VEC)|evallg(ru+1); return roo;
}

GEN
get_mul_table(GEN x,GEN bas,GEN *ptT)
{
  long i,j,k, n = lgef(x)-3;
  GEN T,sym,mul=cgetg(n*n+1,t_MAT);
  
  for (j=1; j<=n*n; j++) mul[j]=lgetg(n+1,t_COL);
  if (ptT)
  {
    sym = polsym(x,n-1); T=cgetg(n+1,t_MAT); *ptT = T;
    for (j=1; j<=n; j++) T[j]=lgetg(n+1,t_COL);
  }
  for (i=1; i<=n; i++)
    for (j=i; j<=n; j++)
    {
      GEN t = gres(gmul((GEN)bas[j],(GEN)bas[i]), x);
      GEN a = (GEN)mul[j+(i-1)*n];
      GEN b = (GEN)mul[i+(j-1)*n];
      long l = lgef(t)-1;
      for (k=1; k<l ; k++) a[k] = b[k] = t[k+1];
      for (   ; k<=n; k++) a[k] = b[k] = zero;
      if (ptT) coeff(T,i,j) = coeff(T,j,i) = (long)quicktrace(t,sym);
    }
  return mul;
}

/* Initialize the number field defined by the polynomial x (in variable v)
 * flag & nf_REGULAR
 *    regular behaviour (no different).
 * flag & nf_DIFFERENT
 *    compute the different.
 * flag & nf_SMALL
 *    compute only nf[1] (pol), nf[2] (signature), nf[5][3] (T2) and
 *    nf[7] (integer basis), the other components are filled with gzero.
 * flag & nf_REDUCE
 *    try a polred first.
 * flag & nf_PARTIAL
 *    do a partial polred, not a polredabs
 * flag & nf_ORIG
 *    do a polred and return [nfinit(x),Mod(a,red)], where
 *    Mod(a,red)=Mod(v,x) (i.e return the base change).
 */

/* here x can be a polynomial, an nf or a bnf */
GEN
initalgall0(GEN x, long flag, long prec)
{
  GEN lead = NULL,nf,ro,bas,mul,mat,M,MC,T,p1,rev,dK,dx,index,fa,res;
  long av=avma,n,i,r1,r2,ru,PRECREG;

  if (DEBUGLEVEL) timer2();
  if (typ(x)==t_POL)
  {
    n=lgef(x)-3; if (n<=0) err(constpoler,"initalgall0");
    check_pol_int(x);
    if (gisirreducible(x) == gzero) err(redpoler,"nfinit");
    if (!gcmp1((GEN)x[n+2]))
    {
      x = pol_to_monic(x,&lead);
      if (!(flag & nf_SMALL))
      {
        if (!(flag & nf_REDUCE))
          err(warner,"non-monic polynomial. Result of the form [nf,c].");
        flag = flag | nf_REDUCE | nf_ORIG;
      }
    }
    bas=allbase4(x,0,&dK,&fa);
    if (DEBUGLEVEL) msgtimer("round4");
    dx = discsr(x); r1 = sturm(x);
  }
  else
  {
    i = lg(x);
    if (typ(x) == t_VEC && i<=4 && i>=3 && typ(x[1])==t_POL)
    { /* polynomial + integer basis */
      bas=(GEN)x[2]; x=(GEN)x[1]; n=lgef(x)-3;
      if (gisirreducible(x) == gzero) err(redpoler,"nfinit");
      dx=discsr(x);
      if (typ(bas) == t_MAT) { mat = bas; bas = mat_to_vecpol(mat,varn(x)); }
      else mat = vecpol_to_mat(bas,n);
      dK = gmul(dx, gsqr(det2(mat)));
      r1 = sturm(x);
    }
    else
    {
      nf=checknf(x);
      bas=(GEN)nf[7]; x=(GEN)nf[1]; n=lgef(x)-3;
      dK=(GEN)nf[3]; dx=mulii(dK, sqri((GEN)nf[4]));
      r1 = itos(gmael(nf,2,1));
    }
    bas[1]=lpolun[varn(x)]; /* it may be gun => SEGV later */
    if (flag & nf_DIFFERENT) fa = factor(absi(dK));
  }
  r2=(n-r1)>>1; ru=r1+r2;
  PRECREG = prec + (expi(dK)>>(TWOPOTBITS_IN_LONG+1))
                 + (long)((sqrt((double)n)+3) / sizeof(long) * 4);
  if (flag & nf_REDUCE)
  {
    nfinit_reduce(flag, &x, &dx, &rev, &bas, r1==n? 0: prec);
    if (DEBUGLEVEL) msgtimer("polred");
  }
  if (!carrecomplet(divii(dx,dK),&index))
    err(bugparier,"nfinit (incorrect discriminant)");

  if (!(flag & nf_SMALL))
  {
    mul = get_mul_table(x,bas, &T);
    if (DEBUGLEVEL) msgtimer("mult. table");
    p1=vecpol_to_mat(bas,n);
  }

  ro=get_roots(x,r1,ru,PRECREG);
  if (DEBUGLEVEL) msgtimer("roots");

  if (flag & nf_ORIG)
  {
    if (!(flag & nf_REDUCE)) err(talker,"bad flag in initalgall0");
    res = cgetg(3,t_VEC);
  }
  nf=cgetg(10,t_VEC);
  nf[1]=lcopy(x);
  nf[2]=lgetg(3,t_VEC);
  mael(nf,2,1)=lstoi(r1);
  mael(nf,2,2)=lstoi(r2);
  nf[3]=lcopy(dK);
  nf[4]=lcopy(index); mat = cgetg(8,t_VEC);
  nf[5]=(long)mat;
  nf[6]=lcopy(ro);
  nf[7]=lcopy(bas);
  if (flag & nf_SMALL) nf[8]=nf[9]=zero;
  else
  {
    nf[8]=linvmat(p1);
    nf[9]=lmul((GEN)nf[8],mul);
  }

  M = make_M(n,ru,bas,ro);
  MC = make_MC(n,r1,ru,M);
  mat[1]=(long)M;
  mat[5]=zero; /* dummy for the different */
  mat[3]=(long)mulmat_real(MC,M);
  if (flag & nf_SMALL)
    mat[2]=mat[4]=mat[6]=mat[7]=zero;
  else
  {
    long av2, av3;
    GEN a2;

    mat[2]=(long)MC;
    mat[4]=lcopy(T);

    av2=avma; p1=ginv(T); av3=avma;
    mat[6] = lpile(av2,av3, gmul(p1,dK));

    av2=avma; p1=content((GEN)mat[6]); a2=gdiv((GEN)mat[6],p1);
    a2 = make_TI(nf,a2,p1); av3=avma;
    /* Ideal basis for discriminant * (inverse of different) */
    mat[7] = lpile(av2,av3,gcopy(a2));
  }
  if (DEBUGLEVEL>=2) msgtimer("matrices");

  if (flag & nf_DIFFERENT)
  {
    mat[5]=(long)differente(nf,fa);
    if (DEBUGLEVEL) msgtimer("different");
  }
  if (!(flag & nf_ORIG)) res = nf;
  else
  {
    res[1]=(long)nf;
    res[2]=lead? ldiv(rev,lead): lcopy(rev);
  }
  return gerepileupto(av,res);
}

GEN
initalgred(GEN x, long prec)
{
  return initalgall0(x,nf_REDUCE|nf_DIFFERENT,prec);
}

GEN
initalgred2(GEN x, long prec)
{
  return initalgall0(x,nf_REDUCE|nf_DIFFERENT|nf_ORIG,prec);
}

GEN
nfinit0(GEN x, long flag,long prec)
{
  switch(flag)
  {
    case 0: return initalgall0(x,nf_DIFFERENT,prec);
    case 1: return initalgall0(x,nf_REGULAR,prec);
    case 2: return initalgall0(x,nf_REDUCE|nf_DIFFERENT,prec);
    case 3: return initalgall0(x,nf_REDUCE|nf_ORIG|nf_DIFFERENT,prec);
    case 4: return initalgall0(x,nf_REDUCE|nf_PARTIAL|nf_DIFFERENT,prec);
    case 5: return initalgall0(x,nf_REDUCE|nf_ORIG|nf_PARTIAL|nf_DIFFERENT,prec);
    case 6: return initalgall0(x,nf_SMALL,prec);
    default: err(flagerr,"nfinit");
  }
  return NULL; /* not reached */
}

GEN
initalg(GEN x, long prec)
{
  return initalgall0(x,nf_DIFFERENT,prec);
}

GEN
nfnewprec(GEN nf, long prec)
{
  long av=avma,i,r1,r2,ru,n,nf_small,tetpil;
  GEN y,pol,ro,bas,MC,mat,M;

  if (typ(nf) != t_VEC) err(talker,"incorrect nf in nfnewprec");
  if (lg(nf) == 11) return bnfnewprec(nf,prec);
  if (prec <= 0) 
  {
    ro = (GEN)nf[6];
    prec = (typ(ro)==t_VEC)?precision((GEN)ro[1]): DEFAULTPREC;
    return (GEN)prec;
  }

  y=cgetg(10,t_VEC);
  for (i=1; i<=4; i++) y[i]=nf[i];
  for (i=6; i<=9; i++) y[i]=nf[i];
  nf_small = gcmp0((GEN)nf[8]);
  pol=(GEN)nf[1]; n=degree(pol);
  r1=itos(gmael(nf,2,1)); r2=itos(gmael(nf,2,2)); ru=r1+r2;
  mat=cgetg(8,t_VEC); y[5]=(long)mat;
  bas=(GEN)nf[7]; ro=get_roots(pol,r1,ru,prec);
  M = make_M(n,ru,bas,ro);
  MC = make_MC(n,r1,ru,M);
  if (nf_small) mat[2]=mat[4]=mat[5]=mat[6]=mat[7]=zero;
  else
  {
    GEN matrices=(GEN)nf[5];
    y[6]=(long)ro;
    mat[2]=(long)MC;
    for (i=4; i<=7; i++) mat[i]=matrices[i];
  }
  mat[1]=(long)M;
  mat[3]=lreal(gmul(MC,M));
  tetpil=avma; return gerepile(av,tetpil,gcopy(y));
}

static long
nf_pm1(GEN y)
{
  long i,l;

  if (!is_pm1(y[1])) return 0;
  l = lg(y);
  for (i=2; i<l; i++)
    if (signe(y[i])) return 0;
  return signe(y[1]);

}

static GEN
is_primitive_root(GEN nf, GEN fa, GEN x, long w)
{
  GEN y, exp = stoi(2), pp = (GEN)fa[1];
  long i,p, l = lg(pp);

  for (i=1; i<l; i++)
  {
    p = itos((GEN)pp[i]);
    exp[2] = w / p; y = element_pow(nf,x,exp);
    if (nf_pm1(y) > 0) /* y = 1 */
    {
      if (p!=2 || !gcmp1(gcoeff(fa,i,2))) return NULL;
      x = gneg_i(x);
    }
  }
  return x;
}

GEN
rootsof1(GEN nf)
{
  long av,tetpil,N,k,i,ws,prec;
  GEN algun,p1,y,R1,d,list,w;

  y=cgetg(3,t_VEC); av=avma; nf=checknf(nf);
  R1=gmael(nf,2,1); algun=gmael(nf,8,1);
  if (signe(R1))
  {
    y[1]=deux;
    y[2]=lneg(algun); return y;
  }
  N=lgef(nf[1])-3; prec=gprecision((GEN)nf[6]);
#ifdef LONG_IS_32BIT
  if (prec < 10) prec = 10;
#else
  if (prec < 6) prec = 6;
#endif
  for (i=1; ; i++)
  {
    p1 = fincke_pohst(gmael(nf,5,3),stoi(N),stoi(1000),1,prec,NULL);
    if (p1) break;
    if (i == MAXITERPOL) err(accurer,"rootsof1");
    prec=(prec<<1)-2;
    if (DEBUGLEVEL) err(warnprec,"rootsof1",prec);
    nf=nfnewprec(nf,prec);
  }
  if (itos(ground((GEN)p1[2])) != N) err(bugparier,"rootsof1 (bug1)");
  w=(GEN)p1[1]; ws = itos(w);
  if (ws == 2)
  {
    y[1]=deux; avma=av;
    y[2]=lneg(algun); return y;
  }

  d = decomp(w); list = (GEN)p1[3]; k = lg(list);
  for (i=1; i<k; i++)
  {
    p1 = (GEN)list[i];
    p1 = is_primitive_root(nf,d,p1,ws);
    if (p1)
    {
      tetpil=avma;
      y[2]=lpile(av,tetpil,gcopy(p1));
      y[1]=lstoi(ws); return y;
    }
  }
  err(bugparier,"rootsof1");
  return NULL; /* not reached */
}

/*******************************************************************/
/*                                                                 */
/*                     DEDEKIND ZETA FUNCTION                      */
/*                                                                 */
/*******************************************************************/

ulong smulss(ulong x, ulong y, ulong *rem);

static GEN
dirzetak0(GEN nf, long N0)
{
  GEN vect,p1,pol,disc,c,c2;
  long av=avma,i,j,k,lx;
  ulong limk,q,p,rem;
  byteptr d=diffptr;
  long court[] = {evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3),0};

  pol=(GEN)nf[1]; disc=(GEN)nf[4];
  c  = (GEN) gpmalloc((N0+1)*sizeof(long));
  c2 = (GEN) gpmalloc((N0+1)*sizeof(long));
  c2[0]=c[0]=evaltyp(t_VEC) | evallg(N0+1);
  c2[1]=c[1]=1; for (i=2; i<=N0; i++) c[i]=0;
  court[2] = 0;

  while (court[2]<=N0)
  {
    court[2] += *d++; if (! *d) err(primer1);
    if (smodis(disc,court[2])) /* court does not divide index */
      { vect = (GEN) simplefactmod(pol,court)[1]; lx=lg(vect); }
    else
    {
      p1=primedec(nf,court); lx=lg(p1); vect=cgetg(lx,t_COL);
      for (i=1; i<lx; i++) vect[i]=mael(p1,i,4);
    }
    for (j=1; j<lx; j++)
    {
      p1=powgi(court, (GEN)vect[j]); /* p1 = court^f */
      if (cmpis(p1,N0) <= 0)
      {
        q=p=p1[2]; limk=N0/q;
        for (k=2; k<=N0; k++) c2[k]=c[k];
        while (q<=N0)
        {
          for (k=1; k<=limk; k++) c2[k*q] += c[k];
          q = smulss(q,p,&rem);
          if (rem) break;
          limk /= p;
        }
        p1=c; c=c2; c2=p1;
      }
    }
    avma=av;
    if (DEBUGLEVEL>6) fprintferr(" %ld",court[2]);
  }
  if (DEBUGLEVEL>6) { fprintferr("\n"); flusherr(); }
  free(c2); return c;
}

GEN
dirzetak(GEN nf, GEN b)
{
  GEN z,c;
  long i;

  if (typ(b)!=t_INT) err(talker,"not an integer type in dirzetak");
  if (signe(b)<=0) return cgetg(1,t_VEC);
  nf = checknf(nf);
  if (is_bigint(b)) err(talker,"too many terms in dirzetak");
  c = dirzetak0(nf,itos(b));
  i = lg(c); z=cgetg(i,t_VEC);
  for (i-- ; i; i--) z[i]=lstoi(c[i]);
  free(c); return z;
}

GEN
initzeta(GEN pol, long prec)
{
  GEN nfz,nf,alpha,beta,mu,gr1,gr2,gru,p1,p2,cst,A0,c0,c1,c2,eps,coef;
  GEN limx,bnf,resi,zet,C,coeflog,racpi,aij,tabj,colzero, *tabcstn, *tabcstni;
  GEN c_even,ck_even,c_odd,ck_odd,serie_even,serie_odd,serie_exp,Pi;
  long N0,imin,imax,r1,r2,ru,R,N,i,j,k,n, av,av2,tetpil;
  long court[] = {evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3),0};
  stackzone *zone, *zone0, *zone1;

  /*************** Calcul du residu et des constantes ***************/
  eps=gmul2n(gun,-bit_accuracy(prec)-6); p1=dbltor(0.5);
  nfz=cgetg(10,t_VEC);
  bnf=buchinit(pol,p1,p1,prec+1); prec=(prec<<1)-1;
  Pi = mppi(prec); racpi=gsqrt(Pi,prec);

  /* Nb de classes et regulateur */
  nf=(GEN)bnf[7]; N=lgef(nf[1])-3;
  gr1=gmael(nf,2,1); gr2=gmael(nf,2,2);
  r1=itos(gr1); r2=itos(gr2); ru=r1+r2; R=ru+2;
  av=avma; p1=(GEN)bnf[8]; p2 = gmul(gmul2n(gmael(p1,1,1),r1), (GEN)p1[2]);
  tetpil = avma; resi=gerepile(av,tetpil,gdiv(p2, gmael(p1,4,1)));

  /* Calcul de N0 */
  cst = cgetr(prec); av = avma;
  mu = gadd(gmul2n(gr1,-1),gr2);
  alpha = gmul2n(stoi(ru+1),-1);
  beta = gpui(gdeux,gmul2n(gr1,-1),DEFAULTPREC);
  A0 = gmul2n(gpowgs(mu,R),r1);
  A0 = gmul(A0,gpowgs(gmul2n(Pi,1),1-ru));
  A0 = gsqrt(A0,DEFAULTPREC);

  c1 = gmul(mu,gpui(beta,ginv(mu),DEFAULTPREC));
  c0 = gdiv(gmul(A0,gpowgs(gmul2n(Pi,1),ru-1)),mu);
  c0 = gmul(c0,gpui(c1,gneg_i(alpha),DEFAULTPREC));
  c2 = gdiv(alpha,mu);

  p1 = glog(gdiv(c0,eps),DEFAULTPREC);
  limx = gdiv(gsub(glog(p1,DEFAULTPREC),glog(c1,DEFAULTPREC)),
              gadd(c2,gdiv(p1,mu)));
  limx = gmul(gpui(gdiv(c1,p1),mu,DEFAULTPREC),
              gadd(gun,gmul(alpha,limx)));
  p1 = gsqrt(absi((GEN)nf[3]),prec);
  p2 = gmul2n(gpowgs(racpi,N),r2);
  gaffect(gdiv(p1,p2), cst);

  av = avma; p1 = gfloor(gdiv(cst,limx)); N0 = p1[2];
  if (is_bigint(p1) || N0 > 10000000)
    err(talker,"discriminant too large for initzeta, sorry");
  if (DEBUGLEVEL>=2)
    { fprintferr("\ninitzeta:\nN0 = %ld\n",N0); flusherr(); timer2(); }

  /* Calcul de imax */

  imin=1; imax=1400;
  p1 = gmul(gpowgs(gmul2n(racpi,1),r2),gpowgs(stoi(5),r1));
  p1 = gdiv(p1,gmul(gmul(gsqrt(limx,DEFAULTPREC),gmul2n(eps,4)),
                         gpowgs(racpi,3)));
  while (imax-imin >= 4)
  {
    long itest = (imax+imin)>>1;
    p2 = gmul(gpowgs(mpfactr(itest,DEFAULTPREC),r2),gpowgs(limx,itest));
    p2 = gmul(p2,gpowgs(mpfactr(itest/2,DEFAULTPREC),r1));
    if (gcmp(p2,p1) >= 0) imax=itest; else imin=itest;
  }
  imax -= (imax & 1); avma = av;
  if (DEBUGLEVEL>=2) { fprintferr("imax = %ld\n",imax); flusherr(); }
  /* Tableau des i/cst (i=1 a N0) */

  i = prec*N0;
  zone  = switch_stack(NULL,i + 2*(N0+1) + 6*prec);
  zone1 = switch_stack(NULL,2*i);
  zone0 = switch_stack(NULL,2*i);
  switch_stack(zone,1);
  tabcstn  = (GEN*) cgetg(N0+1,t_VEC);
  tabcstni = (GEN*) cgetg(N0+1,t_VEC);
  p1 = ginv(cst);
  for (i=1; i<=N0; i++) { tabcstni[i] = gun; tabcstn[i] = mulsr(i,p1); }
  switch_stack(zone,0);

  /********** Calcul des coefficients a(i,j) independants de s **********/

  zet=cgetg(R,t_VEC); zet[1] = lmpeuler(prec);
  for (i=2; i<R; i++) zet[i] = (long)gzeta(stoi(i),prec);

  aij=cgetg(imax+1,t_VEC);
  for (i=1; i<=imax; i++) aij[i] = lgetg(R,t_VEC);

  c_even = realun(prec);
  c_even = gmul2n(c_even,r1);
  c_odd = gmul(c_even,gpowgs(racpi,r1));
  if (ru&1) c_odd=gneg_i(c_odd);
  ck_even=cgetg(R,t_VEC); ck_odd=cgetg(r2+2,t_VEC);
  for (k=1; k<R; k++)
  {
    ck_even[k]=lmul((GEN)zet[k],gadd(gr2,gmul2n(gr1,-k)));
    if (k&1) ck_even[k]=lneg((GEN)ck_even[k]);
  }
  gru=stoi(ru);
  for (k=1; k<=r2+1; k++)
  {
    ck_odd[k]=lmul((GEN)zet[k], gsub(gru, gmul2n(gr1,-k)));
    if (k&1) ck_odd[k]=lneg((GEN)ck_odd[k]);
    ck_odd[k]=ladd(gru,(GEN)ck_odd[k]);
  }
  ck_odd[1]=lsub((GEN)ck_odd[1],gmul(gr1,glog(gdeux,prec)));
  serie_even =cgetg(ru+3,t_SER); serie_odd=cgetg(r2+3,t_SER);
  serie_even[1] = serie_odd[1] = evalsigne(1)+evalvalp(1);
  i=0;

  while (i<imax/2)
  {
    for (k=1; k<R; k++)
      serie_even[k+1]=ldivgs((GEN)ck_even[k],k);
    serie_exp=gmul(c_even,gexp(serie_even,0));
    p1=(GEN)aij[2*i+1];
    for (j=1; j<R; j++) p1[j]=serie_exp[ru+3-j];

    for (k=1; k<=r2+1; k++)
      serie_odd[k+1]=ldivgs((GEN)ck_odd[k],k);
    serie_exp=gmul(c_odd,gexp(serie_odd,0));
    p1=(GEN)aij[2*i+2];
    for (j=1; j<=r2+1; j++) p1[j]=serie_exp[r2+3-j];
    for (   ; j<R; j++) p1[j]=zero;
    i++;

    c_even = gdiv(c_even,gmul(gpowgs(stoi(i),ru),gpowgs(stoi(2*i-1),r2)));
    c_odd  = gdiv(c_odd, gmul(gpowgs(stoi(i),r2),gpowgs(stoi(2*i+1),ru)));
    c_even = gmul2n(c_even,-r2);
    c_odd  = gmul2n(c_odd,r1-r2);
    if (r1&1) { c_even=gneg_i(c_even); c_odd=gneg_i(c_odd); }
    p1 = gr2; p2 = gru;
    for (k=1; k<R; k++)
    {
      p1=gdivgs(p1,2*i-1); p2=gdivgs(p2,2*i);
      ck_even[k] = ladd((GEN)ck_even[k], gadd(p1,p2));
    }
    p1 = gru; p2 = gr2;
    for (k=1; k<=r2+1; k++)
    {
      p1=gdivgs(p1,2*i+1); p2=gdivgs(p2,2*i);
      ck_odd[k] = ladd((GEN)ck_odd[k], gadd(p1,p2));
    }
  }
  tetpil=avma; aij=gerepile(av,tetpil,gcopy(aij));
  if (DEBUGLEVEL>=2) msgtimer("a(i,j)");
  p1=cgetg(5,t_VEC);
  p1[1]=lstoi(r1); p1[2]=lstoi(r2); p1[3]=lstoi(imax); p1[4]=(long)bnf;
  nfz[1]=(long)p1;
  nfz[2]=(long)resi;
  nfz[5]=(long)cst;
  nfz[6]=llog(cst,prec);
  nfz[7]=(long)aij;

  /************* Calcul du nombre d'ideaux de norme donnee *************/
  coef = dirzetak0(nf,N0); tabj = cgetg(N0+1,t_MAT);
  if (DEBUGLEVEL>=2) msgtimer("coef");
  colzero=cgetg(ru+2,t_COL); for (j=1; j<=ru+1; j++) colzero[j]=zero;
  for (i=1; i<=N0; i++)
    if (coef[i])
    {
      GEN t = cgetg(ru+2,t_COL);
      p1 = glog((GEN)tabcstn[i],prec); setsigne(p1, -signe(p1));
      t[1] = lstoi(coef[i]); t[2] = lmul((GEN)t[1],p1);
      for (j=2; j<=ru; j++)
      {
        long av2 = avma; p2 = gmul((GEN)t[j],p1);
        t[j+1] = (long)gerepileuptoleaf(av2, divrs(p2,j));
      }
      /* tabj[n,j]=coef(n)*ln(c/n)^(j-1)/(j-1)! */
      tabj[i] = (long)t;
    }
    else tabj[i]=(long)colzero;
  if (DEBUGLEVEL>=2) msgtimer("a(n)");

  coeflog=cgetg(N0+1,t_VEC); coeflog[1]=zero;
  for (i=2; i<=N0; i++)
    if (coef[i])
    {
      court[2]=i; p1=glog(court,prec);
      setsigne(p1,-1); coeflog[i]=(long)p1;
    }
    else coeflog[i] = zero;
  if (DEBUGLEVEL>=2) msgtimer("log(n)");

  nfz[3]=(long)tabj;
  p1 = cgetg(N0+1,t_VECSMALL);
  for (i=1; i<=N0; i++) p1[i] = coef[i];
  nfz[8]=(long)p1;
  nfz[9]=(long)coeflog;

  /******************** Calcul des coefficients Cik ********************/

  C = cgetg(ru+1,t_MAT);
  for (k=1; k<=ru; k++) C[k] = lgetg(imax+1,t_COL);
  av2 = avma;
  for (i=1; i<=imax; i++)
  {
    stackzone *z;
    for (k=1; k<=ru; k++)
    {
      p1 = NULL;
      for (n=1; n<=N0; n++)
        if (coef[n])
          for (j=1; j<=ru-k+1; j++)
          {
            p2 = gmul(tabcstni[n],
                      gmul(gmael(aij,i,j+k), gmael(tabj,n,j)));
            p1 = p1? gadd(p1,p2): p2;
          }
      coeff(C,i,k) = p1? (long)gerepileupto(av2,p1): zero;
      av2 = avma;
    }
    /* use a parallel stack */
    z = i&1? zone1: zone0;
    switch_stack(z, 1);
    for (n=1; n<=N0; n++)
      if (coef[n]) tabcstni[n] = mpmul(tabcstni[n],tabcstn[n]);
    /* come back */
    switch_stack(z, 0);
  }
  nfz[4] = (long) C;
  if (DEBUGLEVEL>=2) msgtimer("Cik");
  gunclone(aij);
  free((void*)zone); free((void*)zone1); free((void*)zone0);
  free((void*)coef); return nfz;
}

GEN
gzetakall(GEN nfz, GEN s, long flag, long prec2)
{
  GEN resi,C,cst,cstlog,coeflog,cs,coef;
  GEN lambd,gammas,gammaunmoins,gammas2,gammaunmoins2,var1,var2;
  GEN p1,unmoins,gexpro,gar,val,valm,valk,valkm;
  long ts = typ(s), r1,r2,ru,imax,i,j,k,N0,sl,prec,bigprec, av = avma;

  if (typ(nfz)!=t_VEC || lg(nfz)!=10 || typ(nfz[1]) != t_VEC)
    err(talker,"not a zeta number field in zetakall");
  if (! is_intreal_t(ts) && ts != t_COMPLEX && ! is_frac_t(ts))
    err(typeer,"gzetakall");
  resi=(GEN)nfz[2]; C=(GEN)nfz[4]; cst=(GEN)nfz[5];
  cstlog=(GEN)nfz[6]; coef=(GEN)nfz[8]; coeflog=(GEN)nfz[9];
  r1  =itos(gmael(nfz,1,1));
  r2  =itos(gmael(nfz,1,2));
  imax=itos(gmael(nfz,1,3));
  N0=lg(coef)-1; ru=r1+r2;
  /* from initzeta. Certainly excessive, at least if LONG_IS_64BIT */
  bigprec = min(precision(cst), (prec2<<1) - 1);
  prec = prec2+1;

  if (ts==t_COMPLEX && gcmp0(gimag(s))) { s=greal(s); ts = typ(s); }
  if (ts==t_REAL && !signe(gfrac(s))) { s=mptrunc(s); ts = t_INT; }
  if (ts==t_INT)
  {
    sl=itos(s);
    if (sl==1) err(talker,"s = 1 is a pole (gzetakall)");
    if (sl==0)
    {
      if (flag) err(talker,"s = 0 is a pole (gzetakall)");
      if (ru == 1)
      {
        if (r1) return gneg(ghalf);
        return gneg(resi);
      }
      return gzero;
    }
    if (sl<0 && (r2 || !odd(sl)))
    {
      if (!flag) return gzero;
      s = subsi(1,s); sl = 1-sl;
    }
    unmoins=subsi(1,s);
    lambd = gdiv(resi, mulis(s,sl-1));
    gammas2=ggamma(gmul2n(s,-1),prec);
    gar=gpowgs(gammas2,r1);
    cs=gexp(gmul(cstlog,s),prec); 	
    val=s; valm=unmoins;
    if (sl<0) /* r2 = 0 && odd(sl) */
    {
      gammaunmoins2=ggamma(gmul2n(unmoins,-1),prec);
      var1=var2=gun;
      for (i=2; i<=N0; i++)
	if (coef[i])
	{
          gexpro=gexp(gmul((GEN)coeflog[i],s),bigprec);
	  var1 = gadd(var1,gmulsg(coef[i],gexpro));
	  var2 = gadd(var2,gdivsg(coef[i],gmulsg(i,gexpro)));
	}
      lambd=gadd(lambd,gmul(gmul(var1,cs),gar));
      lambd=gadd(lambd,gmul(gmul(var2,gdiv(cst,cs)),
			    gpowgs(gammaunmoins2,r1)));
      for (i=1; i<=imax; i+=2)
      {
	valk  = val;
        valkm = valm;
	for (k=1; k<=ru; k++)	
	{
          p1 = gcoeff(C,i,k);
	  lambd = gsub(lambd,gdiv(p1,valk )); valk  = mulii(val,valk);
	  lambd = gsub(lambd,gdiv(p1,valkm)); valkm = mulii(valm,valkm);
	}
	val  = addis(val, 2);
        valm = addis(valm,2);
      }
    }
    else
    {
      GEN tabj=(GEN)nfz[3], aij=(GEN)nfz[7];

      gar = gmul(gar,gpowgs(ggamma(s,prec),r2));
      var1=var2=gzero;
      for (i=1; i<=N0; i++)
	if (coef[i])
	{
	  gexpro=gexp(gmul((GEN)coeflog[i],s),bigprec);
	  var1 = gadd(var1,gmulsg(coef[i],gexpro));
          if (sl <= imax)
          {
            p1=gzero;
            for (j=1; j<=ru+1; j++)
              p1 = gadd(p1, gmul(gmael(aij,sl,j), gmael(tabj,i,j)));
            var2=gadd(var2,gdiv(p1,gmulsg(i,gexpro)));
          }
	}
      lambd=gadd(lambd,gmul(gmul(var1,cs),gar));
      lambd=gadd(lambd,gmul(var2,gdiv(cst,cs)));
      for (i=1; i<=imax; i++)
      {
	valk  = val;
        valkm = valm;
        if (i == sl)
          for (k=1; k<=ru; k++)
          {	
            p1 = gcoeff(C,i,k);
            lambd = gsub(lambd,gdiv(p1,valk)); valk = mulii(val,valk);
          }
        else
	for (k=1; k<=ru; k++)
	{	
            p1 = gcoeff(C,i,k);
            lambd = gsub(lambd,gdiv(p1,valk )); valk  = mulii(val,valk); 
            lambd = gsub(lambd,gdiv(p1,valkm)); valkm = mulii(valm,valkm);
	}
	val  = addis(val,1);
        valm = addis(valm,1);
      }
    }
  }
  else
  {
    GEN Pi = mppi(prec);
    if (is_frac_t(ts))
      s = gmul(s, realun(bigprec));
    else
      s = gprec_w(s, bigprec);

    unmoins = gsub(gun,s);
    lambd = gdiv(resi,gmul(s,gsub(s,gun)));
    gammas = ggamma(s,prec);
    gammas2= ggamma(gmul2n(s,-1),prec);
    gar = gmul(gpowgs(gammas,r2),gpowgs(gammas2,r1));
    cs = gexp(gmul(cstlog,s),prec);
    var1 = gmul(Pi,s);
    gammaunmoins = gdiv(Pi,gmul(gsin(var1,prec),gammas));
    gammaunmoins2= gdiv(gmul(gmul(gsqrt(Pi,prec),gpui(gdeux,gsub(s,gun),prec)),
                             gammas2),
                        gmul(gcos(gmul2n(var1,-1),prec),gammas));
    var1 = var2 = gun;
    for (i=2; i<=N0; i++)
      if (coef[i])
      {
        gexpro = gexp(gmul((GEN)coeflog[i],s),bigprec);
	var1 = gadd(var1,gmulsg(coef[i], gexpro));
	var2 = gadd(var2,gdivsg(coef[i], gmulsg(i,gexpro)));
      }
    lambd = gadd(lambd,gmul(gmul(var1,cs),gar));
    lambd = gadd(lambd,gmul(gmul(gmul(var2,gdiv(cst,cs)),
	 		         gpowgs(gammaunmoins,r2)),
                            gpowgs(gammaunmoins2,r1)));
    val  = s;
    valm = unmoins;
    for (i=1; i<=imax; i++)
    {
      valk  = val;
      valkm = valm;
      for (k=1; k<=ru; k++)
      {
        p1 = gcoeff(C,i,k);
	lambd = gsub(lambd,gdiv(p1,valk )); valk  = gmul(val, valk);
	lambd = gsub(lambd,gdiv(p1,valkm)); valkm = gmul(valm,valkm);
      }
      if (r2)
      {
	val  = gadd(val, gun);
        valm = gadd(valm,gun);
	}
      else
      {
	val  = gadd(val, gdeux);
        valm = gadd(valm,gdeux); i++;
      }
    }
  }
  if (!flag) lambd = gdiv(lambd,gmul(gar,cs)); /* zetak */
  return gerepileupto(av, gprec_w(lambd, prec2));
}

GEN
gzetak(GEN nfz, GEN s, long prec)
{
  return gzetakall(nfz,s,0,prec);
}

GEN
glambdak(GEN nfz, GEN s, long prec)
{
  return gzetakall(nfz,s,1,prec);
}