/* $Id: subcyclo.c,v 1.21 2002/07/04 12:25:04 bill Exp $

Copyright (C) 2000  The PARI group.

This file is part of the PARI/GP package.

PARI/GP is free software; you can redistribute it and/or modify it under the
terms of the GNU General Public License as published by the Free Software
Foundation. It is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY WHATSOEVER.

Check the License for details. You should have received a copy of it, along
with the package; see the file 'COPYING'. If not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */

#include "pari.h"

extern GEN vandermondeinversemod(GEN L, GEN T, GEN den, GEN mod);
extern GEN supnorm(GEN L, long prec);

/*************************************************************************/
/**                                                                     **/
/**              Routines for handling subgroups of (Z/nZ)^*            **/
/**              without requiring discrete logarithms.                 **/
/**                                                                     **/
/*************************************************************************/

/* Subgroups are [gen,ord,bits] where 
 * gen is a vecsmall of generators 
 * ord is theirs relative orders
 * bits is a bit vector of the elements, of length(n).
 */

 /*The algorithm is similar to testpermutation*/
void
znstar_partial_coset_func(long n, GEN H, void (*func)(void *data,long c)
    , void *data, long d, long c)
{
  GEN gen = (GEN) H[1];
  GEN ord = (GEN) H[2];
  GEN cache = vecsmall_const(d,c);
  long i, j, card = 1;

  (*func)(data,c);
  for (i = 1; i <= d; i++) card *= ord[i];
  for(i=1; i<card; i++)
  {
    long k, m = i;
    for(j=1; j<d && m%ord[j]==0 ;j++) m /= ord[j];
    cache[j] = mulssmod(cache[j],gen[j],n);
    for (k=1; k<j; k++) cache[k] = cache[j];
    (*func)(data, cache[j]);
  }
}

void
znstar_coset_func(long n, GEN H, void (*func)(void *data,long c)
    , void *data, long c)
{
  znstar_partial_coset_func(n, H, func,data, lg(H[1])-1, c);
}

/*Add the element of the bitvec of the coset c modulo the subgroup of H
 * generated by the first d generators to the bitvec bits.*/

void
znstar_partial_coset_bits_inplace(long n, GEN H, GEN bits, long d, long c)
{
  gpmem_t ltop=avma;
  znstar_partial_coset_func(n,H, (void (*)(void *,long)) &bitvec_set,
      (void *) bits, d, c);
  avma=ltop;
}

void
znstar_coset_bits_inplace(long n, GEN H, GEN bits, long c)
{
  znstar_partial_coset_bits_inplace(n, H, bits, lg(H[1])-1, c);
}

GEN
znstar_partial_coset_bits(long n, GEN H, long d, long c)
{
  GEN bits=bitvec_alloc(n);
  znstar_partial_coset_bits_inplace(n,H,bits,d,c);
  return bits;
}
  
/*Compute the bitvec of the elements of the  subgroup of H generated by the
 * first d generators.
 */

GEN
znstar_coset_bits(long n, GEN H, long c)
{
  return znstar_partial_coset_bits(n, H, lg(H[1])-1, c);
}

/*Compute the bitvec of the elements of the  subgroup of H generated by the
 * first d generators.*/

GEN
znstar_partial_bits(long n, GEN H, long d)
{
  return znstar_partial_coset_bits(n, H, d, 1);
}
/*Compute the bitvec of the elements of H.*/

GEN
znstar_bits(long n, GEN H)
{
  return znstar_partial_bits(n,H,lg(H[1])-1);
}

/*Compute the subgroup of (Z/nZ)^* generated by the elements of
 * the vecsmall V.
 */

GEN
znstar_generate(long n, GEN V)
{
  gpmem_t ltop=avma;
  GEN res=cgetg(4,t_VEC);
  GEN gen=cgetg(lg(V),t_VECSMALL);
  GEN ord=cgetg(lg(V),t_VECSMALL);
  GEN bits;
  long i,r=0;
  res[1]=(long)gen;
  res[2]=(long)ord;
  bits=znstar_partial_bits(n,res,r);
  for(i=1;i<lg(V);i++)
  {
    long v=V[i];
    long g=v;
    long o=0;
    while(!bitvec_test(bits,g))
    {
      g=mulssmod(g,v,n);
      o++;
    }
    if (o)
    {
      r++;
      gen[r]=v;
      ord[r]=o+1;
      cgiv(bits); 
      bits=znstar_partial_bits(n,res,r);
    }
  }
  setlg(gen,r+1);
  setlg(ord,r+1);
  res[3]=(long)bits;
  return gerepilecopy(ltop,res);
}

/* Return the lists of element of H.
 * This can be implemented with znstar_coset_func instead.
 */

GEN
znstar_elts(long n, GEN H)
{
  long card=group_order(H);
  GEN gen=(GEN)H[1], ord=(GEN)H[2];
  GEN sg = cgetg(1 + card, t_VECSMALL);
  long k, j, l;
  sg[1] = 1;
  for (j = 1, l = 1; j < lg(gen); j++)
  {
    int     c = l * (ord[j] - 1);
    for (k = 1; k <= c; k++)	/* I like it */
      sg[++l] = mulssmod(sg[k], gen[j], n);
  }
  vecsmall_sort(sg);
  return sg;
}

/* Take a znstar H and n dividing the modulus of H.
 * Output H reduced to modulus n */

GEN
znstar_reduce_modulus(GEN H, long n)
{
  gpmem_t ltop=avma;
  GEN gen=cgetg(lg(H[1]),t_VECSMALL);
  long i;
  for(i=1; i < lg(gen); i++)
    gen[i] = mael(H,1,i)%n;
  return gerepileupto(ltop, znstar_generate(n,gen));
}
/*Compute conductor of H*/

long znstar_conductor(long n, GEN H)
{
  gpmem_t ltop=avma;
  int i,j;
  GEN F;
  long cnd=n;
  F = decomp_small(n); 
  for(i=lg((GEN)F[1])-1;i>0;i--)
  {
    long p=coeff(F,i,1);
    long e=coeff(F,i,2);
    long q=n;
    if (DEBUGLEVEL>=4)
      fprintferr("SubCyclo: testing %ld^%ld\n",p,e);
    for (  ; e>=1; e--)
    {
      long z = 1;
      q /= p;
      for (j = 1; j < p; j++)
      {
	z += q;
	if (!bitvec_test((GEN) H[3],z) && cgcd(z,n)==1)
          break;
      } 
      if ( j < p )
      {
	if (DEBUGLEVEL>=4)
	  fprintferr("SubCyclo: %ld not found\n",z);
	break;
      }
      cnd /= p;
      if (DEBUGLEVEL>=4)
	fprintferr("SubCyclo: new conductor:%ld\n",cnd);
    }
  } 
  if (DEBUGLEVEL>=6)
    fprintferr("SubCyclo: conductor:%ld\n",cnd);
  avma=ltop;
  return cnd;
}

/* Calcule les orbites d'un sous-groupe de Z/nZ donne par un
 * generateur ou d'un ensemble de generateur donne par un vecteur. 
 */
GEN
znstar_cosets(long n, long phi_n, GEN H)
{
  long    k;
  long    c = 0;
  long    card   = group_order(H);
  long    index  = phi_n/card;
  GEN     cosets = cgetg(index+1,t_VECSMALL);
  gpmem_t ltop = avma;
  GEN     bits   = bitvec_alloc(n);
  for (k = 1; k <= index; k++)
  {
    for (c++ ; bitvec_test(bits,c) || cgcd(c,n)!=1; c++);
    cosets[k]=c;
    znstar_coset_bits_inplace(n, H, bits, c);
  }
  avma=ltop;
  return cosets;
}


/*************************************************************************/
/**                                                                     **/
/**                     znstar/HNF interface                            **/
/**                                                                     **/
/*************************************************************************/

/* Convert a true znstar output by znstar to a `small znstar'
 */

GEN
znstar_small(GEN zn)
{
  GEN Z=cgetg(4,t_VEC);
  Z[1]=licopy(gmael3(zn,3,1,1));
  Z[2]=(long) gtovecsmall((GEN)zn[2]);
  Z[3]=(long) lift((GEN)zn[3]);
  return Z;
}


/* Compute generators for the subgroup of (Z/nZ)* given in HNF. 
 */
GEN
znstar_hnf_generators(GEN Z, GEN M)
{
  long l = lg(M);
  GEN gen=cgetg(l, t_VECSMALL);
  gpmem_t ltop=avma;
  GEN zgen= (GEN) Z[3];
  long n = itos((GEN) Z[1]);
  GEN m = stoi(n);
  long j,h;
  for (j = 1; j < l; j++)
  {
    gen[j] = 1;
    for (h = 1; h < l; h++)
      gen[j] = mulssmod(gen[j], 
          itos(powmodulo((GEN) zgen[h], gmael(M,j,h),m)),n);
  }
  avma=ltop;
  return gen;
}

GEN
znstar_hnf(GEN Z, GEN M)
{
  return znstar_generate(itos((GEN)Z[1]),znstar_hnf_generators(Z,M));
}

GEN
znstar_hnf_elts(GEN Z, GEN H)
{
  gpmem_t ltop=avma;
  GEN G=znstar_hnf(Z,H);
  long n=itos((GEN)Z[1]);
  GEN list=znstar_elts(n,G);
  return gerepileupto(ltop,list);
}

/*************************************************************************/
/**                                                                     **/
/**                     subcyclo                                        **/
/**                                                                     **/
/*************************************************************************/

static GEN gscycloconductor(GEN g, long n, long flag)
{
  if (flag==2)
  {
    GEN V=cgetg(3,t_VEC);
    V[1]=lcopy(g);
    V[2]=lstoi(n);
    return V;
  }
  return g;
}

static long 
lift_check_modulus(GEN H, long n)
{
  long t=typ(H);
  long h;
  switch(t)
  {
    case t_INTMOD:
      if (cmpsi(n,(GEN)H[1]))
	err(talker,"wrong modulus in galoissubcyclo");
      H = (GEN)H[2];
    case t_INT:
      h=smodis(H,n);
      if (cgcd(h,n)!=1)
	err(talker,"generators must be prime to conductor in galoissubcyclo");
      return h;
  }
  err(talker,"wrong type in galoissubcyclo");
  return 0;/*not reached*/
}

GEN subcyclo_complex_bound(gpmem_t ltop, GEN V, long prec)
{
  GEN pol = roots_to_pol(V,0);
  GEN vec = gtovec(greal(pol));
  GEN borne = ceil_safe(supnorm(vec,prec));
  return gerepileupto(ltop,borne);
}

GEN subcyclo_complex_cyclic(long n, long d, long m ,long z, long g, GEN powz, long prec)
{
  GEN V=cgetg(d+1,t_VEC);
  long base=1;
  long i,k;
  for (i=1;i<=d;i++,base=mulssmod(base,z,n))
  {
    gpmem_t ltop=avma;
    long ex=base;
    GEN s=gzero;
    (void)new_chunk(2*prec + 3);
    for (k=0; k<m; k++, ex = mulssmod(ex,g,n))
      s=gadd(s,(GEN)powz[ex]);
    avma=ltop;
    V[i]=lcopy(s);
  }
  return V;
}

/* Newton sums mod le. if le==NULL, works with complex instead */
GEN
subcyclo_cyclic(long n, long d, long m ,long z, long g, GEN powz, GEN le)
{
  GEN V=cgetg(d+1,t_VEC);
  long base=1;
  long i,k;
  long lle=le?lg(le)*2+1:2*lg(powz[1])+3;/*Assume dvmdii use lx+ly space*/
  for (i=1;i<=d;i++,base=mulssmod(base,z,n))
  {
    gpmem_t ltop=avma;
    long ex=base;
    GEN s=gzero;
    (void)new_chunk(lle); /* HACK */
    for (k=0; k<m; k++, ex = mulssmod(ex,g,n))
      s=gadd(s,(GEN)powz[ex]);
    avma=ltop;
    V[i]=le?lmodii(s,le):lcopy(s);
  }
  return V;
}

struct _subcyclo_orbits_s
{
  GEN powz;
  GEN *s;
  ulong count;
  gpmem_t ltop;
};

void
_subcyclo_orbits(struct _subcyclo_orbits_s *data, long k)
{
  GEN powz = data->powz;
  GEN *s = data->s;
  
  if (!data->count) data->ltop= avma;
  *s = gadd(*s,(GEN)powz[k]);
  data->count++;
  if ((data->count & 0xffUL) == 0)
    *s = gerepileupto(data->ltop, *s);
}

/* Newton sums mod le. if le==NULL, works with complex instead */
GEN
subcyclo_orbits(long n, GEN H, GEN O, GEN powz, GEN le)
{
  long i, d=lg(O);
  GEN V=cgetg(d,t_VEC);
  struct _subcyclo_orbits_s data;
  long lle=le?lg(le)*2+1:2*lg(powz[1])+3;/*Assume dvmdii use lx+ly space*/
  data.powz = powz;
  for(i=1; i<d; i++)
  {
    GEN s = gzero;
    gpmem_t av = avma;
    (void)new_chunk(lle);
    data.count = 0;
    data.s     = &s;
    znstar_coset_func(n, H, (void (*)(void *,long)) _subcyclo_orbits,
      (void *) &data, O[i]);
    avma = av; /* HACK */
    V[i] = le?lmodii(s,le):lcopy(s);
  }
  return V;
}

GEN 
subcyclo_start(long n, long d, long o, GEN borne, long *ptr_val,long *ptr_l)
{
  gpmem_t av;
  GEN l,le,z;
  long i;
  long e,val;
  if (DEBUGLEVEL >= 1) (void)timer2();
  l=stoi(n+1);e=1;
  while(!isprime(l)) 
  { 
    l=addis(l,n);
    e++;
  }
  if (DEBUGLEVEL >= 4)
    fprintferr("Subcyclo: prime l=%Z\n",l);
  av=avma;
  if (!borne)
  {
    /*Borne utilise': 
      Vecmax(Vec((x+o)^d)=max{binome(d,i)*o^i ;1<=i<=d} 
     */
    i=d-(1+d)/(1+o);
    borne=mulii(binome(stoi(d),i),gpowgs(stoi(o),i));
  }
  if (DEBUGLEVEL >= 4)
    fprintferr("Subcyclo: borne=%Z\n",borne);
  val=logint(shifti(borne,2),l,NULL);
  avma=av;
  if (DEBUGLEVEL >= 4)
    fprintferr("Subcyclo: val=%ld\n",val);
  le=gpowgs(l,val);
  z=lift(gpowgs(gener(l),e));
  z=padicsqrtnlift(gun,stoi(n),z,l,val);
  if (DEBUGLEVEL >= 1)
    msgtimer("padicsqrtnlift.");
  *ptr_val=val;
  *ptr_l=itos(l);
  return gmodulcp(z,le);
}

GEN
subcyclo_complex_roots(long n, long real, long prec)
{
  GEN powz, z = exp_Ir(divrs(Pi2n(1, prec), n)); /* = e_n(1) */
  long i, k = (n+3)>>1;

  powz = cgetg(n,t_VEC);
  powz[1] = (long)z;
  for (i=2; i<k; i++) powz[i] = lmul(z,(GEN)powz[i-1]);
  if (real) /* 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]);
  return powz;
}

GEN
subcyclo_roots(long n, GEN zl)
{
  GEN le=(GEN) zl[1];
  GEN z=(GEN) zl[2];
  long lle=lg(le)*3; /*Assume dvmdii use lx+ly space*/
  long i;
  GEN powz = cgetg(n,t_VEC);
  powz[1] = (long) z;
  for (i=2; i<n; i++)
  {
    gpmem_t av=avma;
    GEN p1;
    (void)new_chunk(lle); /* HACK */
    p1 = mulii(z,(GEN)powz[i-1]);
    avma=av;
    powz[i] = lmodii(p1,le);
  }
  return powz;
}

GEN
galoiscyclo(long n, long v)
{
  ulong ltop=avma;
  GEN grp,G;
  GEN z, le;
  long val,l;
  GEN L;
  long i,j,k;
  GEN zn=znstar(stoi(n));
  long card=itos((GEN) zn[1]);
  GEN gen=lift((GEN)zn[3]);
  GEN ord=gtovecsmall((GEN)zn[2]);
  GEN elts;
  z=subcyclo_start(n,card/2,2,NULL,&val,&l);
  le=(GEN) z[1];
  z=(GEN) z[2];
  L = cgetg(1+card,t_VEC);
  L[1] = (long) z;
  for (j = 1, i = 1; j < lg(gen); j++)
  {
    int     c = i * (ord[j] - 1);
    for (k = 1; k <= c; k++)	/* I like it */
      L[++i] = (long) powmodulo((GEN)L[k],(GEN)gen[j],le);
  }
  G=abelian_group(ord);
  elts = group_elts(G, card); /*not stack clean*/
  grp = cgetg(9, t_VEC);
  grp[1] = (long) cyclo(n,v);
  grp[2] = lgetg(4,t_VEC); 
  mael(grp,2,1) = lstoi(l);
  mael(grp,2,2) = lstoi(val);
  mael(grp,2,3) = licopy(le);
  grp[3] = lcopy(L);
  grp[4] = (long) vandermondeinversemod(L, (GEN) grp[1], gun, le);
  grp[5] = un;
  grp[6] = lcopy(elts);
  grp[7] = lcopy((GEN)G[1]);
  grp[8] = lcopy((GEN)G[2]);
  return gerepileupto(ltop,grp);
}

/* Convert a bnrinit(Q,n) to a znstar(n)
 * complex is set to 0 if the bnr is real and to 1 if it is complex.
 * Not stack clean 
 */
GEN bnrtozn(GEN bnr, long *complex)
{
  GEN zk;
  GEN gen;
  GEN cond;
  long l2;
  long i;
  GEN p3;         /* vec */
  GEN res;
  checkbnrgen(bnr);
  zk = (GEN) bnr[5];
  gen = (GEN) zk[3];
  /*cond is the finite part of the conductor
   * complex is the infinite part*/
  cond = gcoeff(gmael3(bnr,2,1,1), 1, 1);
  *complex = signe(gmael4(bnr,2,1,2,1));
  l2 = lg(gen);
  res= cgetg(4,t_VEC); 
  res[1]=zk[1];
  res[2]=zk[2];
  p3 = cgetg(l2, t_VEC);
  for (i = 1; i < l2; ++i)
  {
    GEN x=(GEN) gen[i];
    if (typ(x) == t_MAT)
      x = gcoeff(x, 1, 1);
    else if (typ(x) == t_COL)
      x = (GEN) x[1];
    p3[i] = (long) gmodulcp(mpabs(x), cond);
  }
  res[3] = (long) p3;
  return res;
}

GEN 
galoissubcyclo(GEN N, GEN sg, long flag, long v)
{
  gpmem_t ltop=avma,av;
  GEN H, V;
  long i;
  GEN O;
  GEN Z=NULL;
  GEN B,zl,L,T,le,powz;
  long val,l;
  long n, cnd, complex=1;
  long card, phi_n;
  if (flag<0 || flag>2) err(flagerr,"galoissubcyclo");
  if ( v==-1 ) v=0;
  if (!sg) sg=gun;
  switch(typ(N))
  {
    case t_INT:
      n=itos(N);
      if ( n<1 ) err(arither2);
      break;
    case t_VEC:
      if (lg(N)==7)
        N=bnrtozn(N,&complex);
      if (lg(N)==4)
      {
        Z=N;
        if (typ(Z[3])!=t_VEC)
          err(typeer,"galoissubcyclo");
        if (lg(Z[3])==1)
          n=1;
        else
        {
          if (typ(gmael(Z,3,1))!= t_INTMOD)
#ifdef NETHACK_MESSAGES
            err(talker,"You have transgressed!");
#else
            err(talker,"Please do not try to break PARI with ridiculously counterfeit data. Thanks!");
#endif
          n=itos(gmael3(Z,3,1,1));
        }
        break;
      }
    default: /*fall through*/
      err(typeer,"galoissubcyclo");
      return NULL;/*Not reached*/
  }
  if (n==1) {avma=ltop; return polx[v];}
  switch(typ(sg))
  {
     case t_INTMOD: case t_INT: 
      V=cgetg(2,t_VECSMALL); 
      V[1]=lift_check_modulus(sg,n);
      break;
    case t_VECSMALL:
      V=gcopy(sg);
      for (i=1;i<lg(V);i++)
        if (V[i]<0)
          V[i]=mulssmod(-V[i],n-1,n);
      break;
    case t_VEC:
    case t_COL:
      V=cgetg(lg(sg),t_VECSMALL);
      for(i=1;i<lg(sg);i++)
        V[i] = (long)lift_check_modulus((GEN)sg[i],n);
      break;
    case t_MAT:/*Fall through*/
      {
        if (lg(sg) == 1 || lg(sg) != lg(sg[1]))
          err(talker,"not a HNF matrix in galoissubcyclo");
        if (!Z)
          err(talker,"N must be a bnrinit or a znstar if H is a matrix in galoissubcyclo");
        if ( lg(Z[2]) != lg(sg) || lg(Z[3]) != lg(sg))
          err(talker,"Matrix of wrong dimensions in galoissubcyclo");
        V = znstar_hnf_generators(znstar_small(Z),sg);
      }
      break;
    default:
      err(typeer,"galoissubcyclo");
      return NULL;/*Not reached*/
  }
  if (!complex) /*Add complex conjugation*/
    V=vecsmall_append(V,n-1);
  H = znstar_generate(n,V);
  /* compute the complex/real status
   * it is real iff z -> conj(z)=z^-1=z^(n-1) is in H
   */
  if (DEBUGLEVEL >= 6)
  {
    fprintferr("Subcyclo: elements:");
    for (i=1;i<n;i++)
      if (bitvec_test((GEN)H[3],i))
        fprintferr(" %ld",i);
    fprintferr("\n");
  }
  complex = !bitvec_test((GEN) H[3],n-1);
  if (DEBUGLEVEL >= 6)
    fprintferr("Subcyclo: complex=%ld\n",complex);
  if (DEBUGLEVEL >= 1) (void)timer2();
  cnd = znstar_conductor(n,H);
  if (DEBUGLEVEL >= 1)
    msgtimer("znstar_conductor");
  if ( flag == 1 )  { avma=ltop; return stoi(cnd); }
  if (n != cnd)
  {
    H = znstar_reduce_modulus(H, cnd);
    n = cnd;
  }
  card = group_order(H);
  phi_n= itos(phi(stoi(n)));
  if ( card==phi_n )
  {
    avma=ltop;
    if (flag==3) return galoiscyclo(n,v);
    return gscycloconductor(cyclo(n,v),n,flag); 
  }
  O = znstar_cosets(n, phi_n, H);
  if (DEBUGLEVEL >= 1)
    msgtimer("znstar_cosets");
  if (DEBUGLEVEL >= 6)
    fprintferr("Subcyclo: orbits=%Z\n",O);
  if (DEBUGLEVEL >= 4)
    fprintferr("Subcyclo: %ld orbits with %ld elements each\n",phi_n/card,card);
  av=avma;
  powz=subcyclo_complex_roots(n,!complex,3);
  L=subcyclo_orbits(n,H,O,powz,NULL);
  B=subcyclo_complex_bound(av,L,3);
  zl=subcyclo_start(n,phi_n/card,card,B,&val,&l);
  powz=subcyclo_roots(n,zl);
  le=(GEN) zl[1];
  L=subcyclo_orbits(n,H,O,powz,le);
  T=FpV_roots_to_pol(L,le,v);
  T=FpX_center(T,le);
  return gerepileupto(ltop,gscycloconductor(T,n,flag));
}

GEN
subcyclo(long n, long d, long v)
{
  gpmem_t ltop=avma;
  long o,p,al,r,g,gd;
  GEN fa,G;
  GEN zl,L,T,le;
  long l,val;
  GEN B,powz;
  if (v<0) v = 0;
  if (d==1) return polx[v];
  if (d<=0 || n<=0) err(typeer,"subcyclo");
  if ((n & 3) == 2) n >>= 1;
  if (n == 1 || d >= n) err(talker,"degree does not divide phi(n) in subcyclo");
  fa = decomp(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(talker,"non-cyclic case in polsubcyclo: use galoissubcyclo instead");
  avma=ltop;
  r = cgcd(d,n); /* = p^(v_p(d))*/
  n = r*p;
  o = n-r; /* = phi(n) */
  if (o == d) return cyclo(n,v);
  if (o % d) err(talker,"degree does not divide phi(n) in subcyclo");
  o /= d;
  if (p==2)
  {
    GEN pol = powgi(polx[v],gdeux); pol[2]=un; /* replace gzero */
    return pol; /* = x^2 + 1 */
  }
  G=gener(stoi(n));
  g=itos((GEN)G[2]);
  gd=itos((GEN)gpowgs(G,d)[2]);
  avma=ltop;
  powz=subcyclo_complex_roots(n,(o&1)==0,3);
  L=subcyclo_cyclic(n,d,o,g,gd,powz,NULL);
  B=subcyclo_complex_bound(ltop,L,3);
  zl=subcyclo_start(n,d,o,B,&val,&l);
  le=(GEN)zl[1];
  powz=subcyclo_roots(n,zl);
  if (DEBUGLEVEL >= 6)
    msgtimer("subcyclo_roots"); 
  L=subcyclo_cyclic(n,d,o,g,gd,powz,le);
  if (DEBUGLEVEL >= 6)
    msgtimer("subcyclo_cyclic"); 
  T=FpV_roots_to_pol(L,le,v);
  if (DEBUGLEVEL >= 6)
    msgtimer("roots_to_pol"); 
  T=FpX_center(T,le);
  return gerepileupto(ltop,T);
}

GEN polsubcyclo(long n, long d, long v)
{
  gpmem_t ltop=avma;
  GEN L, Z=znstar(stoi(n));
  /*subcyclo is twice faster but Z must be cyclic*/
  if (lg(Z[2]) == 2 && divise((GEN)Z[1],stoi(d)))
  {
    avma=ltop; 
    return subcyclo(n, d, v);
  }
  L=subgrouplist((GEN) Z[2], _vec(stoi(d)));
  if (lg(L) == 2)
    return gerepileupto(ltop, galoissubcyclo(Z, (GEN) L[1], 0, v));
  else
  {
    GEN V=cgetg(lg(L),t_VEC);
    long i;
    for (i=1; i< lg(V); i++)
      V[i] = (long) galoissubcyclo(Z, (GEN) L[i], 0, v);
    return gerepileupto(ltop,V);
  }
}