OpenXM_contrib/pari-2.2/src/modules/stark.c - diff

Return to stark.c CVS log

Up to [local] / OpenXM_contrib / pari-2.2 / src / modules

Diff for /OpenXM_contrib/pari-2.2/src/modules/Attic/stark.c between version 1.1 and 1.2

-version 1.1, 2001/10/02 11:17:12
+version 1.2, 2002/09/11 07:27:06
 Line 15  Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 Line 15  Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 Line 15  Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  /*******************************************************************/
  /*                                                                 */
- /*                  COMPUTATION OF STARK UNITS                     */
+ /*        COMPUTATION OF STARK UNITS OF TOTALLY REAL FIELDS        */
- /*                    OF TOTALLY REAL FIELDS                       */
  /*                                                                 */
  /*******************************************************************/
  #include "pari.h"
-Line 26  Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+Line 25  Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 Line 26  Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 Line 25  Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  #define ADD_PREC   (DEFAULTPREC-2)*3
  extern GEN roots_to_pol_intern(GEN L, GEN a, long v, int plus);
+ extern GEN bnrGetSurj(GEN bnr1, GEN bnr2);
+ /* ComputeCoeff */
+ typedef struct {
+   GEN L0, L1, L11, L2; /* VECSMALL of p */
+   GEN *L1ray, *L11ray; /* precomputed isprincipalray(pr), pr | p */
+   GEN *rayZ; /* precomputed isprincipalray(i), i < condZ */
+   long condZ; /* generates cond(bnr) \cap Z, assumed small */
+ } LISTray;
+ /* Char evaluation */
+ typedef struct {
+   long ord;
+   GEN *val, chi;
+ } CHI_t;
+ /* RecCoeff */
+ typedef struct {
+   GEN M, beta, B, U, nB;
+   long v, G, N;
+ } RC_data;
  /********************************************************************/
  /*                    Miscellaneous functions                       */
  /********************************************************************/
+ /* exp(2iPi/den), assume den a t_INT */
- /* Compute the image of logelt by chi as a complex number if flag = 0,
+ GEN
-    otherwise as a polmod, see InitChar in part 3 */
+ InitRU(GEN den, long prec)
+ {
+   GEN c,s, z;
+   if (egalii(den, gdeux)) return stoi(-1);
+   gsincos(divri(gmul2n(mppi(prec),1), den), &s, &c, prec);
+   z = cgetg(3, t_COMPLEX);
+   z[1] = (long)c;
+   z[2] = (long)s; return z;
+ }
+ /* Compute the image of logelt by chi as a complex number
+    see InitChar in part 3 */
  static GEN
- ComputeImagebyChar(GEN chi, GEN logelt, long flag)
+ ComputeImagebyChar(GEN chi, GEN logelt)
  {
-   GEN gn = gmul((GEN)chi[1], logelt), x = (GEN)chi[flag? 4: 2];
+   GEN gn = gmul((GEN)chi[1], logelt), x = (GEN)chi[2];
    long d = itos((GEN)chi[3]), n = smodis(gn, d);
    /* x^d = 1 and, if d even, x^(d/2) = -1 */
    if ((d & 1) == 0)
-Line 47  ComputeImagebyChar(GEN chi, GEN logelt, long flag)
+Line 77  ComputeImagebyChar(GEN chi, GEN logelt, long flag)
 Line 47  ComputeImagebyChar(GEN chi, GEN logelt, long flag)
 Line 77  ComputeImagebyChar(GEN chi, GEN logelt, long flag)
    return gpowgs(x, n);
  }
+ /* return n such that C(elt) = z^n */
+ static long
+ EvalChar_n(CHI_t *C, GEN logelt)
+ {
+   GEN n = gmul(C->chi, logelt);
+   return smodis(n, C->ord);
+ }
+ /* return C(elt) */
+ static GEN
+ EvalChar(CHI_t *C, GEN logelt)
+ {
+   return C->val[EvalChar_n(C, logelt)];
+ }
+ static void
+ init_CHI(CHI_t *c, GEN CHI, GEN z)
+ {
+   long i, d = itos((GEN)CHI[3]);
+   GEN *v = (GEN*)new_chunk(d);
+   v[1] = z;
+   for (i=2; i<d; i++)
+     v[i] = gmul(v[i-1], z);
+   v[0] = gmul(v[i-1], z);
+   c->chi = (GEN)CHI[1];
+   c->ord = d;
+   c->val = v;
+ }
+ /* as t_COMPLEX */
+ static void
+ init_CHI_alg(CHI_t *c, GEN CHI) {
+   GEN z = gmodulcp(polx[0], cyclo(itos((GEN)CHI[3]), 0));
+   init_CHI(c,CHI,z);
+ }
+ /* as t_POLMOD */
+ static void
+ init_CHI_C(CHI_t *c, GEN CHI) {
+   init_CHI(c,CHI, (GEN)CHI[2]);
+ }
  /* Compute the conjugate character */
  static GEN
  ConjChar(GEN chi, GEN cyc)
-Line 59  ConjChar(GEN chi, GEN cyc)
+Line 129  ConjChar(GEN chi, GEN cyc)
 Line 59  ConjChar(GEN chi, GEN cyc)
 Line 129  ConjChar(GEN chi, GEN cyc)
        p1[i] = zero;
      else
        p1[i] = lsubii((GEN)cyc[i], (GEN)chi[i]);
    return p1;
  }
- /* compute the next element for FindEltofGroup */
+ typedef struct {
- static GEN
+   long r; /* rank = lg(gen) */
- NextEltofGroup(GEN cyc, long l, long adec)
+   GEN j; /* current elt is gen[1]^j[1] ... gen[r]^j[r] */
- {
+   GEN cyc; /* t_VECSMALL of elementary divisors */
-   GEN p1;
+ } GROUP_t;
-   long dj, j;
-   p1 = cgetg(l + 1, t_COL);
+ static int
+ NextElt(GROUP_t *G)
-   for (j = 1; j <= l; j++)
+ {
+   long i = 1;
+   if (G->r == 0) return 0; /* no more elt */
+   while (++G->j[i] == G->cyc[i]) /* from 0 to cyc[i]-1 */
    {
-     dj = itos((GEN)cyc[j]);
+     G->j[i] = 0;
-     p1[j] = lstoi(adec%dj);
+     if (++i > G->r) return 0; /* no more elt */
-     adec /= dj;
    }
+   return i; /* we have multiplied by gen[i] */
-   return p1;
  }
  /* Compute all the elements of a group given by its SNF */
  static GEN
- FindEltofGroup(long order, GEN cyc)
+ EltsOfGroup(long order, GEN cyc)
  {
-   long l, i;
+   long i;
    GEN rep;
+   GROUP_t G;
-   l = lg(cyc)-1;
+   G.cyc = gtovecsmall(cyc);
+   G.r = lg(cyc)-1;
+   G.j = vecsmall_const(G.r, 0);
    rep = cgetg(order + 1, t_VEC);
+   rep[order] = (long)small_to_col(G.j);
-   for  (i = 1; i <= order; i++)
+   for  (i = 1; i < order; i++)
-     rep[i] = (long)NextEltofGroup(cyc, l, i);
+   {
+     (void)NextElt(&G);
+     rep[i] = (long)small_to_col(G.j);
+   }
    return rep;
  }
-Line 104  FindEltofGroup(long order, GEN cyc)
+Line 180  FindEltofGroup(long order, GEN cyc)
 Line 104  FindEltofGroup(long order, GEN cyc)
 Line 180  FindEltofGroup(long order, GEN cyc)
  static GEN
  ComputeLift(GEN dataC)
  {
-   long order, i, av = avma;
+   long order, i;
+   gpmem_t av = avma;
    GEN cyc, surj, eltq, elt;
    order = itos((GEN)dataC[1]);
    cyc   = (GEN)dataC[2];
    surj  = (GEN)dataC[3];
-   eltq = FindEltofGroup(order, cyc);
+   eltq = EltsOfGroup(order, cyc);
    elt = cgetg(order + 1, t_VEC);
-Line 121  ComputeLift(GEN dataC)
+Line 198  ComputeLift(GEN dataC)
 Line 121  ComputeLift(GEN dataC)
 Line 198  ComputeLift(GEN dataC)
    return gerepileupto(av, elt);
  }
- /* Let bnr1, bnr2 be such that mod(bnr2) | mod(bnr1), compute the
-    matrix of the surjective map Cl(bnr1) ->> Cl(bnr2) */
  static GEN
- GetSurjMat(GEN bnr1, GEN bnr2)
+ get_chic(GEN chi, GEN cyc)
  {
-   long nbg, i;
+   long i, l = lg(chi);
-   GEN gen, M;
+   GEN chic = cgetg(l, t_VEC);
+   for (i = 1; i < l; i++)
-   gen = gmael(bnr1, 5, 3);
+     chic[i] = ldiv((GEN)chi[i], (GEN)cyc[i]);
-   nbg = lg(gen) - 1;
+   return chic;
-   M = cgetg(nbg + 1, t_MAT);
-   for (i = 1; i <= nbg; i++)
-     M[i] = (long)isprincipalray(bnr2, (GEN)gen[i]);
-   return M;
  }
  /* A character is given by a vector [(c_i), z, d, pm] such that
     chi(id) = z ^ sum(c_i * a_i) where
       a_i= log(id) on the generators of bnr
-      z  = exp(2i * Pi / d)
+      z  = exp(2i * Pi / d) */
-      pm = z as a polmod */
  static GEN
- get_Char(GEN chi, long prec)
+ get_Char(GEN chic, long prec)
  {
-   GEN p2, d, _2ipi = gmul(gi, shiftr(mppi(prec), 1));
+   GEN d, C;
-   p2 = cgetg(5, t_VEC); d = denom(chi);
+   C = cgetg(4, t_VEC);
-   p2[1] = lmul(d, chi);
+   d = denom(chic);
-   if (egalii(d, gdeux))
+   C[1] = lmul(d, chic);
-     p2[2] = lstoi(-1);
+   C[2] = (long)InitRU(d, prec);
-   else
+   C[3] = (long)d; return C;
-     p2[2] = lexp(gdiv(_2ipi, d), prec);
-   p2[3] = (long)d;
-   p2[4] = lmodulcp(polx[0], cyclo(itos(d), 0));
-   return p2;
  }
- /* Let chi a character defined over bnr and primitif over bnrc,
+ /* Let chi a character defined over bnr and primitive over bnrc,
     compute the corresponding primitive character and the vectors of
-    prime ideals dividing bnr but not bnr. Returns NULL if bnr = bnrc */
+    prime ideals dividing bnr but not bnrc. Returns NULL if bnr = bnrc */
  static GEN
  GetPrimChar(GEN chi, GEN bnr, GEN bnrc, long prec)
  {
-   long nbg, i, j, l, av = avma, nd;
+   long nbg, i, j, l, nd;
-   GEN gen, cyc, U, chic, M, s, p1, cond, condc, p2, nf;
+   gpmem_t av = avma;
+   GEN s, chic, cyc, U, M, p1, cond, condc, p2, nf;
    GEN prdiff, Mrc;
-   cond  = gmael(bnr, 2, 1);
+   cond  = gmael(bnr,  2, 1);
    condc = gmael(bnrc, 2, 1);
    if (gegal(cond, condc)) return NULL;
-   gen   = gmael(bnr, 5, 3);
+   cyc   = gmael(bnr, 5, 2); nbg = lg(cyc)-1;
-   nbg   = lg(gen) - 1;
-   cyc   = gmael(bnr, 5, 2);
    Mrc   = diagonal(gmael(bnrc, 5, 2));
    nf    = gmael(bnr, 1, 7);
-   cond  = (GEN)cond[1];
+   M = bnrGetSurj(bnr, bnrc);
-   condc = (GEN)condc[1];
+   U = (GEN)hnfall(concatsp(M, Mrc))[2];
+   l = lg((GEN)M[1]);
-   M  = GetSurjMat(bnr, bnrc);
-   l  = lg((GEN)M[1]);
-   p1 = hnfall(concatsp(M, Mrc));
-   U  = (GEN)p1[2];
    chic = cgetg(l, t_VEC);
    for (i = 1; i < l; i++)
    {
-Line 194  GetPrimChar(GEN chi, GEN bnr, GEN bnrc, long prec)
+Line 252  GetPrimChar(GEN chi, GEN bnr, GEN bnrc, long prec)
 Line 194  GetPrimChar(GEN chi, GEN bnr, GEN bnrc, long prec)
 Line 252  GetPrimChar(GEN chi, GEN bnr, GEN bnrc, long prec)
      for (j = 1; j <= nbg; j++)
      {
        p2 = gdiv((GEN)p1[j], (GEN)cyc[j]);
-       s  = gadd(s,gmul(p2,(GEN)chi[j]));
+       s  = gadd(s, gmul(p2,(GEN)chi[j]));
      }
      chic[i] = (long)s;
    }
+   cond  = (GEN)cond[1];
+   condc = (GEN)condc[1];
    p2 = (GEN)idealfactor(nf, cond)[1];
    l  = lg(p2);
-Line 214  GetPrimChar(GEN chi, GEN bnr, GEN bnrc, long prec)
+Line 274  GetPrimChar(GEN chi, GEN bnr, GEN bnrc, long prec)
 Line 214  GetPrimChar(GEN chi, GEN bnr, GEN bnrc, long prec)
 Line 274  GetPrimChar(GEN chi, GEN bnr, GEN bnrc, long prec)
    return gerepileupto(av,p1);
  }
- /* Let dataCR be a list of characters, compute the image of logelt */
  static GEN
- chiideal(GEN dataCR, GEN logelt, long flag)
- {
-   long j, l = lg(dataCR);
-   GEN rep = cgetg(l, t_VEC);
-   for (j = 1; j < l; j++)
-     rep[j] = (long)ComputeImagebyChar(gmael(dataCR, j, 5), logelt, flag);
-   return rep;
- }
- static GEN
  GetDeg(GEN dataCR)
  {
    long i, l = lg(dataCR);
    GEN degs = cgetg(l, t_VECSMALL);
    for (i = 1; i < l; i++)
-     degs[i] = degpol(gmael4(dataCR, i, 5, 4, 1));
+     degs[i] = itos(phi(gmael3(dataCR, i, 5, 3)));
    return degs;
  }
-Line 244  GetDeg(GEN dataCR)
+Line 291  GetDeg(GEN dataCR)
 Line 244  GetDeg(GEN dataCR)
 Line 291  GetDeg(GEN dataCR)
  static GEN AllStark(GEN data,  GEN nf,  long flag,  long prec);
  static GEN InitChar0(GEN dataD, long prec);
+ static GEN QuadGetST(GEN dataCR, GEN vChar, long prec);
  /* Let A be a finite abelian group given by its relation and let C
     define a subgroup of A, compute the order of A / C, its structure and
-Line 251  static GEN InitChar0(GEN dataD, long prec);
+Line 299  static GEN InitChar0(GEN dataD, long prec);
 Line 251  static GEN InitChar0(GEN dataD, long prec);
 Line 299  static GEN InitChar0(GEN dataD, long prec);
  static GEN
  InitQuotient0(GEN DA, GEN C)
  {
-   long i, l;
+   GEN D, MQ, MrC, H, U, rep;
-   GEN MQ, MrC, H, snf, snf2, rep, p1;
-   l = lg(DA)-1;
    H = gcmp0(C)? DA: C;
    MrC  = gauss(H, DA);
-   snf  = smith2(hnf(MrC));
+   (void)smithall(hnf(MrC), &U, NULL);
-   MQ   = concatsp(gmul(H, (GEN)snf[1]), DA);
+   MQ   = concatsp(gmul(H, U), DA);
-   snf2 = smith2(hnf(MQ));
+   D = smithall(hnf(MQ), &U, NULL);
    rep = cgetg(5, t_VEC);
+   rep[1] = (long)dethnf_i(D);
-   p1  = cgetg(l + 1, t_VEC);
+   rep[2] = (long)mattodiagonal(D);
-   for (i = 1; i <= l; i++)
+   rep[3] = lcopy(U);
-     p1[i] = lcopy(gcoeff((GEN)snf2[3], i, i));
-   rep[1] = (long)dethnf((GEN)snf2[3]);
-   rep[2] = (long)p1;
-   rep[3] = lcopy((GEN)snf2[1]);
    rep[4] = lcopy(C);
    return rep;
-Line 286  static GEN
+Line 327  static GEN
 Line 286  static GEN
 Line 327  static GEN
  InitQuotient(GEN bnr, GEN C)
  {
    GEN Mrm, dataquo = cgetg(3, t_VEC);
-   long av;
+   gpmem_t av;
    dataquo[1] = lcopy(bnr);
    av = avma;  Mrm = diagonal(gmael(bnr, 5, 2));
-Line 296  InitQuotient(GEN bnr, GEN C)
+Line 337  InitQuotient(GEN bnr, GEN C)
 Line 296  InitQuotient(GEN bnr, GEN C)
 Line 337  InitQuotient(GEN bnr, GEN C)
  }
  /* Let s: A -> B given by P, and let DA, DB be resp. the matrix of the
-    relations of A and B, let nbA, nbB be resp. the rank of A and B,
+    relations of A and B, compute the kernel of s. If DA = 0 then A is free */
-    compute the kernel of s. If DA = 0 then A is free */
  static GEN
- ComputeKernel0(GEN P, GEN DA, GEN DB, long nbA, long nbB)
+ ComputeKernel0(GEN P, GEN DA, GEN DB)
  {
-   long rk, av = avma;
+   gpmem_t av = avma;
-   GEN herm, mask1, mask2, U;
+   long nbA = lg(DA)-1, rk;
+   GEN herm, U;
    herm  = hnfall(concatsp(P, DB));
-   rk = nbA + nbB + 1;
+   rk = nbA + lg(DB) - lg(herm[1]);
-   rk -= lg((GEN)herm[1]); /* two steps: bug in pgcc 1.1.3 inlining (IS) */
+   U = (GEN)herm[2];
+   U = vecextract_i(U, 1,rk);
+   U = rowextract_i(U, 1,nbA);
-   mask1 = subis(shifti(gun, nbA), 1);
-   mask2 = subis(shifti(gun, rk), 1);
-   U = matextract((GEN)herm[2], mask1, mask2);
    if (!gcmp0(DA)) U = concatsp(U, DA);
    return gerepileupto(av, hnf(U));
  }
  /* Let m and n be two moduli such that n|m and let C be a congruence
     group modulo n, compute the corresponding congruence group modulo m
-    ie then kernel of the map Clk(m) ->> Clk(n)/C */
+    ie the kernel of the map Clk(m) ->> Clk(n)/C */
  static GEN
  ComputeKernel(GEN bnrm, GEN dataC)
  {
-   long av = avma, i, nbm, nbq;
+   long i, nbm;
+   gpmem_t av = avma;
    GEN bnrn, Mrm, genm, Mrq, mgq, P;
    bnrn = (GEN)dataC[1];
    Mrm  = diagonal(gmael(bnrm, 5, 2));
+   Mrq  = diagonal(gmael(dataC, 2, 2));
    genm = gmael(bnrm, 5, 3);
    nbm  = lg(genm) - 1;
-   Mrq  = diagonal(gmael(dataC, 2, 2));
    mgq  = gmael(dataC, 2, 3);
-   nbq  = lg(mgq) - 1;
    P = cgetg(nbm + 1, t_MAT);
    for (i = 1; i <= nbm; i++)
      P[i] = lmul(mgq, isprincipalray(bnrn, (GEN)genm[i]));
-   return gerepileupto(av, ComputeKernel0(P, Mrm, Mrq, nbm, nbq));
+   return gerepileupto(av, ComputeKernel0(P, Mrm, Mrq));
  }
  /* Let C a congruence group in bnr, compute its subgroups of index 2 as
-Line 346  ComputeKernel(GEN bnrm, GEN dataC)
+Line 384  ComputeKernel(GEN bnrm, GEN dataC)
 Line 346  ComputeKernel(GEN bnrm, GEN dataC)
 Line 384  ComputeKernel(GEN bnrm, GEN dataC)
  static GEN
  ComputeIndex2Subgroup(GEN bnr, GEN C)
  {
-   long nb, i, l, av = avma;
+   gpmem_t av = avma;
-   GEN snf, Mr, U, CU, subgrp, rep, p1;
+   long nb, i;
+   GEN D, Mr, U, T, subgrp;
    disable_dbg(0);
    Mr = diagonal(gmael(bnr, 5, 2));
-   snf = smith2(gmul(ginv(C), Mr));
+   D = smithall(gauss(C, Mr), &U, NULL);
+   T = gmul(C,ginv(U));
+   subgrp  = subgrouplist(D, gdeux);
+   nb = lg(subgrp) - 1;
+   setlg(subgrp, nb); /* skip Id which comes last */
+   for (i = 1; i < nb; i++)
+     subgrp[i] = (long)hnf(concatsp(gmul(T, (GEN)subgrp[i]), Mr));
-   U = ginv((GEN)snf[1]);
+   disable_dbg(-1);
-   l = lg((GEN)snf[3]);
+   return gerepilecopy(av, subgrp);
+ }
-   p1 = cgetg(l, t_VEC);
+ GEN
+ Order(GEN cyc, GEN x)
+ {
+   gpmem_t av = avma;
+   long i, l = lg(cyc);
+   GEN c,o,f = gun;
    for (i = 1; i < l; i++)
-     p1[i] = mael3(snf, 3, i, i);
+   {
+     o = (GEN)cyc[i];
-   subgrp  = subgrouplist(p1, 2);
+     c = mppgcd(o, (GEN)x[i]);
-   nb = lg(subgrp) - 1; CU = gmul(C,U);
+     if (!is_pm1(c)) o = diviiexact(o,c);
+     f = mpppcm(f, o);
-   rep = cgetg(nb, t_VEC);
+   }
-   for (i = 1; i < nb; i++) /* skip Id which comes last */
+   return gerepileuptoint(av, f);
-     rep[i] = (long)hnf(concatsp(gmul(CU, (GEN)subgrp[i]), Mr));
-   disable_dbg(-1);
-   return gerepilecopy(av, rep);
  }
  /* Let pr be a prime (pr may divide mod(bnr)), compute the indexes
-Line 377  ComputeIndex2Subgroup(GEN bnr, GEN C)
+Line 424  ComputeIndex2Subgroup(GEN bnr, GEN C)
 Line 377  ComputeIndex2Subgroup(GEN bnr, GEN C)
 Line 424  ComputeIndex2Subgroup(GEN bnr, GEN C)
  static GEN
  GetIndex(GEN pr, GEN bnr, GEN subgroup)
  {
-   long av = avma, v, lg, i;
+   long v, e, f;
-   GEN bnf, mod, mod0, mpr0, mpr, bnrpr, subpr, M, e, f, dtQ, p1;
+   gpmem_t av = avma;
+   GEN bnf, mod, mod0, bnrpr, subpr, M, dtQ, p1;
    GEN rep, cycpr, cycQ;
    bnf  = (GEN)bnr[1];
    mod  = gmael(bnr, 2, 1);
    mod0 = (GEN)mod[1];
-   /* Compute the part of mod coprime to pr */
    v = idealval(bnf, mod0, pr);
-   mpr0 = idealdivexact(bnf, mod0, idealpow(bnf, pr, stoi(v)));
+   if (v == 0)
-   mpr = cgetg(3, t_VEC);
-   mpr[1] = (long)mpr0;
-   mpr[2] = mod[2];
-   if (gegal(mpr0, mod0))
    {
      bnrpr = bnr;
      subpr = subgroup;
+     e = 1;
    }
    else
    {
+     GEN mpr = cgetg(3, t_VEC);
+     GEN mpr0 = idealdivpowprime(bnf, mod0, pr, stoi(v));
+     mpr[1] = (long)mpr0; /* part of mod coprime to pr */
+     mpr[2] = mod[2];
      bnrpr = buchrayinitgen(bnf, mpr);
      cycpr = gmael(bnrpr, 5, 2);
-     M = GetSurjMat(bnr, bnrpr);
+     M = bnrGetSurj(bnr, bnrpr);
      M = gmul(M, subgroup);
      subpr = hnf(concatsp(M, diagonal(cycpr)));
+     /* e = #(bnr/subgroup) / #(bnrpr/subpr) */
+     e = itos( divii(dethnf_i(subgroup), dethnf_i(subpr)) );
    }
-   /* e = #(bnr/subgroup) / #(bnrpr/subpr) */
-   e = gdiv(det(subgroup), det(subpr));
    /* f = order of [pr] in bnrpr/subpr */
    dtQ  = InitQuotient(bnrpr, subpr);
    p1   = isprincipalray(bnrpr, pr);
    p1   = gmul(gmael(dtQ, 2, 3), p1);
    cycQ = gmael(dtQ, 2, 2);
-   lg = lg(cycQ) - 1;
+   f  = itos( Order(cycQ, p1) );
-   f  = gun;
-   for (i = 1; i <= lg; i++)
-     f = glcm(f, gdiv((GEN)cycQ[i], ggcd((GEN)cycQ[i], (GEN)p1[i])));
-   rep = cgetg(3, t_VEC);
+   rep = cgetg(3, t_VECSMALL);
-   rep[1] = lcopy(e);
+   rep[1] = (long)e;
-   rep[2] = lcopy(f);
+   rep[2] = (long)f;
    return gerepileupto(av, rep);
  }
-Line 432  GetIndex(GEN pr, GEN bnr, GEN subgroup)
+Line 475  GetIndex(GEN pr, GEN bnr, GEN subgroup)
 Line 432  GetIndex(GEN pr, GEN bnr, GEN subgroup)
 Line 475  GetIndex(GEN pr, GEN bnr, GEN subgroup)
  static GEN
  CplxModulus(GEN data, long *newprec, long prec)
  {
-   long av = avma, pr, dprec;
+   long pr, dprec;
+   gpmem_t av = avma, av2;
    GEN nf, cpl, pol, p1;
    nf = gmael3(data, 1, 1, 7);
-Line 450  CplxModulus(GEN data, long *newprec, long prec)
+Line 494  CplxModulus(GEN data, long *newprec, long prec)
 Line 450  CplxModulus(GEN data, long *newprec, long prec)
 Line 494  CplxModulus(GEN data, long *newprec, long prec)
    dprec = DEFAULTPREC;
+   av2 = avma;
    for (;;)
    {
      p1[5] = (long)InitChar0((GEN)data[3], dprec);
-     pol   = AllStark(p1, nf, -1, dprec);
+     pol   = gerepileupto(av2, AllStark(p1, nf, -1, dprec));
      if (!gcmp0(leading_term(pol)))
      {
-       cpl   = mpabs(poldisc0(pol, 0));
+       cpl = gnorml2(gtovec(pol));
        if (!gcmp0(cpl)) break;
      }
      pr = gexpo(pol)>>(TWOPOTBITS_IN_LONG+1);
-Line 488  FindModulus(GEN dataC, long fl, long *newprec, long pr
+Line 533  FindModulus(GEN dataC, long fl, long *newprec, long pr
 Line 488  FindModulus(GEN dataC, long fl, long *newprec, long pr
 Line 533  FindModulus(GEN dataC, long fl, long *newprec, long pr
  {
    long n, i, narch, nbp, maxnorm, minnorm, N, nbidnn, s, c, j, nbcand;
    long limnorm, first = 1, pr;
-   ulong av = avma, av1, av0;
+   gpmem_t av = avma, av1, av0;
    GEN bnr, rep, bnf, nf, f, arch, m, listid, idnormn, bnrm, ImC;
    GEN candD, D, bpr, indpr, sgp, p1, p2, rb;
-Line 503  FindModulus(GEN dataC, long fl, long *newprec, long pr
+Line 548  FindModulus(GEN dataC, long fl, long *newprec, long pr
 Line 503  FindModulus(GEN dataC, long fl, long *newprec, long pr
 Line 548  FindModulus(GEN dataC, long fl, long *newprec, long pr
    for (i = 1; i <= 5; i++) rep[i] = zero;
    /* if cpl < rb, it is not necessary to try another modulus */
-   rb = powgi(gmul(gmael(bnf, 7, 3), det(f)), gmul2n(gmael(bnr, 5, 1), 3));
+   rb = powgi(gmul((GEN)nf[3], det(f)), gmul2n(gmael(bnr, 5, 1), 3));
    bpr = (GEN)idealfactor(nf, f)[1];
    nbp = lg(bpr) - 1;
-   indpr = cgetg(nbp + 1,t_VEC);
+   indpr = cgetg(nbp + 1,t_VECSMALL);
    for (i = 1; i <= nbp; i++)
    {
      p1 = GetIndex((GEN)bpr[i], bnr, sgp);
-     indpr[i] = lmulii((GEN)p1[1], (GEN)p1[2]);
+     indpr[i] = p1[1] * p1[2];
    }
    /* Initialization of the possible infinite part */
-Line 529  FindModulus(GEN dataC, long fl, long *newprec, long pr
+Line 574  FindModulus(GEN dataC, long fl, long *newprec, long pr
 Line 529  FindModulus(GEN dataC, long fl, long *newprec, long pr
 Line 574  FindModulus(GEN dataC, long fl, long *newprec, long pr
       If we cannot find a suitable conductor of norm < limnorm, we stop */
    maxnorm = 50;
    minnorm = 1;
-   limnorm = 200;
+   limnorm = 400;
    if (DEBUGLEVEL >= 2)
      fprintferr("Looking for a modulus of norm: ");
-Line 566  FindModulus(GEN dataC, long fl, long *newprec, long pr
+Line 611  FindModulus(GEN dataC, long fl, long *newprec, long pr
 Line 566  FindModulus(GEN dataC, long fl, long *newprec, long pr
 Line 611  FindModulus(GEN dataC, long fl, long *newprec, long pr
            bnrm = buchrayinitgen(bnf, m);
            p1   = conductor(bnrm, gzero, -1);
            disable_dbg(-1);
-           if (signe(p1))
+           arch[N+1-s] = un;
-           {
+           if (!signe(p1)) continue;
-             /* compute Im(C) in Clk(m)... */
-             ImC = ComputeKernel(bnrm, dataC);
-             /* ... and its subgroups of index 2 */
+           /* compute Im(C) in Clk(m)... */
-             candD  = ComputeIndex2Subgroup(bnrm, ImC);
+           ImC = ComputeKernel(bnrm, dataC);
-             nbcand = lg(candD) - 1;
-             for (c = 1; c <= nbcand; c++)
-             {
-               /* check if m is the conductor */
-               D  = (GEN)candD[c];
-               disable_dbg(0);
-               p1 = conductor(bnrm, D, -1);
-               disable_dbg(-1);
-               if (signe(p1))
-               {
-                 /* check the splitting of primes */
-                 for (j = 1; j <= nbp; j++)
-                 {
-                   p1 = GetIndex((GEN)bpr[j], bnrm, D);
-                   p1 = mulii((GEN)p1[1], (GEN)p1[2]);
-                   if (egalii(p1, (GEN)indpr[j])) break; /* no good */
-                 }
-                 if (j > nbp)
-                 {
-                   p2 = cgetg(6, t_VEC);
-                   p2[1] = lcopy(bnrm);
+           /* ... and its subgroups of index 2 */
-                   p2[2] = lcopy(D);
+           candD  = ComputeIndex2Subgroup(bnrm, ImC);
-                   p2[3] = (long)InitQuotient((GEN)p2[1], (GEN)p2[2]);
+           nbcand = lg(candD) - 1;
-                   p2[4] = (long)InitQuotient((GEN)p2[1], ImC);
+           for (c = 1; c <= nbcand; c++)
+           {
+             /* check if m is the conductor */
+             D  = (GEN)candD[c];
+             disable_dbg(0);
+             p1 = conductor(bnrm, D, -1);
+             disable_dbg(-1);
+             if (!signe(p1)) continue;
-                   p1 = CplxModulus(p2, &pr, prec);
+             /* check the splitting of primes */
+             for (j = 1; j <= nbp; j++)
+             {
+               p1 = GetIndex((GEN)bpr[j], bnrm, D);
+               if (p1[1] * p1[2] == indpr[j]) break; /* no good */
+             }
+             if (j <= nbp) continue;
-                   if (first == 1 || gcmp(p1, (GEN)rep[5]) < 0)
+             p2 = cgetg(6, t_VEC);
-                   {
+             p2[1] = (long)bnrm;
-                     *newprec = pr;
+             p2[2] = (long)D;
-                     for (j = 1; j <= 4; j++) rep[j] = p2[j];
+             p2[3] = (long)InitQuotient(bnrm, D);
-                     rep[5] = (long)p1;
+             p2[4] = (long)InitQuotient(bnrm, ImC);
-                   }
-                   if (!fl || (gcmp(p1, rb) < 0))
+             p1 = CplxModulus(p2, &pr, prec);
-                   {
-                     rep[5] = (long)InitChar0((GEN)rep[3], *newprec);
+             if (first == 1 || gcmp(p1, (GEN)rep[5]) < 0)
-                     return gerepilecopy(av, rep);
+             {
-                   }
+               *newprec = pr;
-                   if (DEBUGLEVEL >= 2)
+               for (j = 1; j <= 4; j++) rep[j] = p2[j];
-                     fprintferr("Trying to find another modulus...");
+               rep[5] = (long)p1;
-                   first--;
+             }
-                 }
-               }
+             if (!fl || gcmp(p1, rb) < 0) goto END; /* OK */
-             }
-           }
+             if (DEBUGLEVEL>1) fprintferr("Trying to find another modulus...");
-           arch[N+1-s] = un;
+             first--;
+           }
          }
          if (first <= bnd)
          {
-           if (DEBUGLEVEL >= 2)
+           if (DEBUGLEVEL>1)
              fprintferr("No, we're done!\nModulus = %Z and subgroup = %Z\n",
                         gmael3(rep, 1, 2, 1), rep[2]);
-           rep[5] = (long)InitChar0((GEN)rep[3], *newprec);
+           goto END; /* OK */
-           return gerepilecopy(av, rep);
          }
        }
      }
-Line 640  FindModulus(GEN dataC, long fl, long *newprec, long pr
+Line 675  FindModulus(GEN dataC, long fl, long *newprec, long pr
 Line 640  FindModulus(GEN dataC, long fl, long *newprec, long pr
 Line 675  FindModulus(GEN dataC, long fl, long *newprec, long pr
      if (maxnorm > limnorm)
        err(talker, "Cannot find a suitable modulus in FindModulus");
    }
+ END:
+   rep[5] = (long)InitChar0((GEN)rep[3], *newprec);
+   return gerepilecopy(av, rep);
  }
  /********************************************************************/
  /*                      2nd part: compute W(X)                      */
  /********************************************************************/
- /* compute W(chi) such that Ld(s,chi) = W(chi) Ld(1 - s, chi*),
+ /* compute the list of W(chi) such that Ld(s,chi) = W(chi) Ld(1 - s, chi*),
-    if flag != 0 do not check the result */
+  * for all chi in LCHI. All chi have the same conductor (= cond(bnr)).
+  * if check == 0 do not check the result */
  static GEN
- ComputeArtinNumber(GEN datachi, long flag, long prec)
+ ComputeArtinNumber(GEN bnr, GEN LCHI, long check, long prec)
  {
-   long av = avma, av2, G, ms, j, i, nz, zcard, q, l, N, lim;
+   long ic, i, j, nz, q, N, nChar = lg(LCHI)-1;
-   GEN chi, nc, dc, p1, cond0, cond1, elts, Msign, umod2, lambda, nf;
+   gpmem_t av = avma, av2, lim;
-   GEN sg, p2, chib, diff, vt, z, idg, mu, idh, zid, zstruc, zgen, zchi;
+   GEN sqrtnc, dc, cond, condZ, cond0, cond1, lambda, nf, T;
-   GEN classe, bnr, beta, s, tr, p3, den, muslambda, Pi, lp1, beta2;
+   GEN cyc, *vN, *vB, diff, vt, idg, mu, idh, zid, *gen, z, *nchi;
+   GEN indW, W, CHI, classe, s0, *s, tr, den, muslambda, beta2, sarch;
+   CHI_t **lC;
+   GROUP_t G;
-   chi   = (GEN)datachi[8];
+   lC = (CHI_t**)new_chunk(nChar + 1);
-   /* trivial case */
+   indW = cgetg(nChar + 1, t_VECSMALL);
-   if (cmpsi(2, (GEN)chi[3]) >= 0) return gun;
+   W = cgetg(nChar + 1, t_VEC);
+   for (ic = 0, i = 1; i <= nChar; i++)
+   {
+     CHI = (GEN)LCHI[i];
+     if (cmpsi(2, (GEN)CHI[3]) >= 0) { W[i] = un; continue; } /* trivial case */
+     ic++; indW[ic] = i;
+     lC[ic] = (CHI_t*)new_chunk(sizeof(CHI_t));
+     init_CHI_C(lC[ic], CHI);
+   }
+   if (!ic) return W;
+   nChar = ic;
-   bnr   = (GEN)datachi[3];
    nf    = gmael(bnr, 1, 7);
    diff  = gmael(nf, 5, 5);
-   cond0 = gmael3(bnr, 2, 1, 1);
+   T     = gmael(nf,5,4);
-   cond1 = gmael3(bnr, 2, 1, 2);
+   cond  =  gmael(bnr, 2, 1);
-   umod2 = gmodulcp(gun, gdeux);
+   cond0 = (GEN)cond[1]; condZ = gcoeff(cond0,1,1);
+   cond1 = (GEN)cond[2];
    N     = degpol(nf[1]);
-   Pi    = mppi(prec);
-   G   = - bit_accuracy(prec) >> 1;
+   sqrtnc  = gsqrt(idealnorm(nf, cond0), prec);
-   nc  = idealnorm(nf, cond0);
    dc  = idealmul(nf, diff, cond0);
    den = idealnorm(nf, dc);
-   z   = gexp(gdiv(gmul2n(gmul(gi, Pi), 1), den), prec);
+   z   = InitRU(den, prec);
    q = 0;
    for (i = 1; i < lg(cond1); i++)
-Line 681  ComputeArtinNumber(GEN datachi, long flag, long prec)
+Line 731  ComputeArtinNumber(GEN datachi, long flag, long prec)
 Line 681  ComputeArtinNumber(GEN datachi, long flag, long prec)
 Line 731  ComputeArtinNumber(GEN datachi, long flag, long prec)
    /* compute a system of elements congru to 1 mod cond0 and giving all
       possible signatures for cond1 */
-   p1 = zarchstar(nf, cond0, cond1, q);
+   sarch = zarchstar(nf, cond0, cond1, q);
-   elts = (GEN)p1[2];
-   Msign = gmul((GEN)p1[3], umod2);
-   ms = lg(elts) - 1;
    /* find lambda in diff.cond such that gcd(lambda.(diff.cond)^-1,cond0) = 1
-      and lambda >(cond1)> 0 */
+      and lambda >> 0 at cond1 */
    lambda = idealappr(nf, dc);
-   sg = zsigne(nf, lambda, cond1);
+   lambda = set_sign_mod_idele(nf, NULL, lambda, cond,sarch);
-   p2 = lift(gmul(Msign, sg));
-   for (j = 1; j <= ms; j++)
-     if (gcmp1((GEN)p2[j])) lambda = element_mul(nf, lambda, (GEN)elts[j]);
    idg = idealdivexact(nf, lambda, dc);
    /* find mu in idg such that idh=(mu) / idg is coprime with cond0 and
-      mu >(cond1)> 0 */
+      mu >> 0 at cond1 */
    if (!gcmp1(gcoeff(idg, 1, 1)))
    {
-     p1 = idealfactor(nf, idg);
+     GEN p1 = idealfactor(nf, idg);
-     p2 = idealfactor(nf, cond0);
+     GEN p2 = idealfactor(nf, cond0);
+     p2[2] = (long)zerocol(lg(p2[1])-1);
+     p1 = concat_factor(p1,p2);
-     l = lg((GEN)p2[1]);
-     for (i = 1; i < l; i++) coeff(p2, i, 2) = zero;
-     p1 = gtrans(concatsp(gtrans(p1), gtrans(p2)));
      mu = idealapprfact(nf, p1);
-     sg = zsigne(nf, mu, cond1);
+     mu = set_sign_mod_idele(nf, NULL,mu, cond,sarch);
-     p2 = lift(gmul(Msign, sg));
-     for (j = 1; j <= ms; j++)
-       if (gcmp1((GEN)p2[j])) mu = element_mul(nf, mu, (GEN)elts[j]);
      idh = idealdivexact(nf, mu, idg);
    }
    else
    {
      mu  = gun;
-     idh = gcopy(idg);
+     idh = idg;
    }
-   muslambda = element_div(nf, mu, lambda);
+   muslambda = gmul(den, element_div(nf, mu, lambda));
    /* compute a system of generators of (Ok/cond)^* cond1-positive */
    zid = zidealstarinitgen(nf, cond0);
+   cyc = gmael(zid, 2, 2);
+   gen = (GEN*)gmael(zid, 2, 3);
+   nz = lg(gen) - 1;
-   zcard  = itos(gmael(zid, 2, 1));
+   nchi = (GEN*)cgetg(nChar+1, t_VEC);
-   zstruc = gmael(zid, 2, 2);
+   for (ic = 1; ic <= nChar; ic++) nchi[ic] = cgetg(nz + 1, t_VECSMALL);
-   zgen   = gmael(zid, 2, 3);
-   nz = lg(zgen) - 1;
-   zchi = cgetg(nz + 1, t_VEC);
    for (i = 1; i <= nz; i++)
    {
-     p1 = (GEN)zgen[i];
+     if (is_bigint(cyc[i]))
-     sg = zsigne(nf, p1, cond1);
+       err(talker,"conductor too large in ComputeArtinNumber");
-     p2 = lift(gmul(Msign, sg));
+     gen[i] = set_sign_mod_idele(nf, NULL,gen[i], cond,sarch);
+     classe = isprincipalray(bnr, gen[i]);
-     for (j = 1; j <= ms; j++)
+     for (ic = 1; ic <= nChar; ic++)
-       if (gcmp1((GEN)p2[j])) p1 = element_mul(nf, p1, (GEN)elts[j]);
+       nchi[ic][i] = (long)EvalChar_n(lC[ic], classe);
-     classe = isprincipalray(bnr, p1);
-     zchi[i] = (long)ComputeImagebyChar(chi, classe, 0);
-     zgen[i] = (long)p1;
    }
    /* Sum chi(beta) * exp(2i * Pi * Tr(beta * mu / lambda) where beta
       runs through the classes of (Ok/cond0)^* and beta cond1-positive */
-   p3 = cgetg(N + 1, t_COL);
+   vt = cgetg(N + 1, t_VEC); /* Tr(w_i) */
-   for (i = 1; i <= N; i++) p3[i] = zero;
+   for (i = 1; i <= N; i++) vt[i] = coeff(T,i,1);
-   vt = cgetg(N + 1, t_VEC);
+   G.cyc = gtovecsmall(cyc);
-   for (i = 1; i <= N; i++)
+   G.r = nz;
-   {
+   G.j = vecsmall_const(nz, 0);
-     p3[i] = un;
-     vt[i] = ltrace(basistoalg(nf, p3));
-     p3[i] = zero;
-   }
-   lp1 = NULL;
+   vN = (GEN*)cgetg(nChar+1, t_VEC);
-   s = gzero;
+   for (ic = 1; ic <= nChar; ic++) vN[ic] = vecsmall_const(nz, 0);
    av2 = avma; lim = stack_lim(av2, 1);
+   vB = (GEN*)cgetg(nz+1, t_VEC);
+   for (i=1; i<=nz; i++) vB[i] = gun;
-   for (i = 0; i < zcard; i++)
+   tr = gmod(gmul(vt, muslambda), den); /* for beta = 1 */
+   s0 = powgi(z, tr);
+   s = (GEN*)cgetg(nChar+1, t_VEC);
+   for (ic = 1; ic <= nChar; ic++) s[ic] = s0;
+   for (;;)
    {
-     p1 = NextEltofGroup(zstruc, nz, i);
+     if (! (i = NextElt(&G)) ) break;
-     /* we test if we can use the previous value
+     vB[i] = FpV_red(element_mul(nf, vB[i], gen[i]), condZ);
-        of beta / chib to compute the next one */
+     for (j=1; j<i; j++) vB[j] = vB[i];
-     /* FIXME: there should be a better way of doing that! */
-     if (!lp1 || !gcmp1(gnorml2(gsub(p1, lp1))))
-     {
-       beta = gun;
-       chib = gun;
-       for (j = 1; j <= nz; j++)
+     for (ic = 1; ic <= nChar; ic++)
-       {
-         if (!gcmp0((GEN)p1[j]))
-         {
-           p2 = element_powmodideal(nf, (GEN)zgen[j], (GEN)p1[j], cond0);
-           beta = element_mulmodideal(nf, beta, p2, cond0);
-           chib = gmul(chib, powgi((GEN)zchi[j], (GEN)p1[j]));
-         }
-       }
-     }
-     else
      {
-       /* we use the fact that NextEltofGroup is, in general,
+       vN[ic][i] = addssmod(vN[ic][i], nchi[ic][i], lC[ic]->ord);
-          obtained by adding 1 to a component of p1 */
+       for (j=1; j<i; j++) vN[ic][j] = vN[ic][i];
-       for (j = 1; j <= nz; j++)
+     }
-         if (!gegal((GEN)p1[j], (GEN)lp1[j]))
-         {
-           beta = element_mulmodideal(nf, beta, (GEN)zgen[j], cond0);
-           chib = gmul(chib, (GEN)zchi[j]);
-         }
-     }
-     lp1 = p1;
+     vB[i]= set_sign_mod_idele(nf, NULL,vB[i], cond,sarch);
-     sg  = zsigne(nf, beta, cond1);
+     beta2 = element_mul(nf, vB[i], muslambda);
-     p2  = lift(gmul(Msign, sg));
-     for (j = 1; j <= ms; j++)
+     tr = gmod(gmul(vt, beta2), den);
-       if (gcmp1((GEN)p2[j]))
+     s0 = powgi(z,tr);
-         beta = element_mul(nf, beta, (GEN)elts[j]);
+     for (ic = 1; ic <= nChar; ic++)
+     {
+       long ind = vN[ic][i];
+       GEN val = (GEN)lC[ic]->val[ind];
+       s[ic] = gadd(s[ic], gmul(val, s0));
+     }
-     beta2 = element_mul(nf, beta, muslambda);
-     tr = gmul(vt, beta2);
-     tr = gmod(gmul(tr, den), den);
-     s = gadd(s, gmul(chib, powgi(z,tr)));
      if (low_stack(lim, stack_lim(av2, 1)))
      {
-       GEN *gptr[5];
+       GEN *gptr[2]; gptr[0] = (GEN*)&s; gptr[1] = (GEN*)&vB;
-       gptr[0] = &s; gptr[1] = &lp1; gptr[2] = &beta; gptr[3] = &chib;
        if (DEBUGMEM > 1) err(warnmem,"ComputeArtinNumber");
-       gerepilemany(av2, gptr, 4);
+       gerepilemany(av2, gptr, 2);
      }
    }
    classe = isprincipalray(bnr, idh);
-   s = gmul(s, ComputeImagebyChar(chi, classe, 0));
+   z = gpowgs(gneg_i(gi),q);
-   s = gdiv(s, gsqrt(nc, prec));
+   for (ic = 1; ic <= nChar; ic++)
+   {
+     s0 = gmul(s[ic], EvalChar(lC[ic], classe));
+     s0 = gdiv(s0, sqrtnc);
+     if (check && - expo(subrs(gnorm(s0), 1)) < bit_accuracy(prec) >> 1)
+       err(bugparier, "ComputeArtinNumber");
+     W[indW[ic]] = lmul(s0, z);
+   }
+   return gerepilecopy(av, W);
+ }
-   p1 = gsubgs(gabs(s, prec), 1);
+ static GEN
+ ComputeAllArtinNumbers(GEN dataCR, GEN vChar, int check, long prec)
+ {
+   const long cl = lg(dataCR) - 1;
+   long ncond = lg(vChar)-1;
+   long jc, k;
+   GEN W = cgetg(cl+1,t_VEC), WbyCond, LCHI;
-   i = expo(p1);
+   for (jc = 1; jc <= ncond; jc++)
-   if ((i > G) && !flag)
+   {
-     err(bugparier, "ComputeArtinNumber");
+     const GEN LChar = (GEN)vChar[jc];
+     const long nChar = lg(LChar)-1;
+     GEN *ldata = (GEN*)vecextract_p(dataCR, LChar);
+     GEN dtcr = ldata[1], bnr = (GEN)dtcr[3];
-   return gerepileupto(av, gmul(s, gpowgs(gneg_i(gi),q)));
+     if (DEBUGLEVEL>1)
+       fprintferr("* Root Number: cond. no %ld/%ld (%ld chars)\n",
+                  jc,ncond,nChar);
+     LCHI = cgetg(nChar + 1, t_VEC);
+     for (k = 1; k <= nChar; k++) LCHI[k] = ldata[k][8];
+     WbyCond = ComputeArtinNumber(bnr, LCHI, check, prec);
+     for (k = 1; k <= nChar; k++) W[LChar[k]] = WbyCond[k];
+   }
+   return W;
  }
  /* compute the constant W of the functional equation of
-Line 842  ComputeArtinNumber(GEN datachi, long flag, long prec)
+Line 878  ComputeArtinNumber(GEN datachi, long flag, long prec)
 Line 842  ComputeArtinNumber(GEN datachi, long flag, long prec)
 Line 878  ComputeArtinNumber(GEN datachi, long flag, long prec)
  GEN
  bnrrootnumber(GEN bnr, GEN chi, long flag, long prec)
  {
-   long av = avma, l, i;
+   long l;
-   GEN cond, condc, bnrc, chic, cyc, d, p1, p2, dtcr, Pi;
+   gpmem_t av = avma;
+   GEN cond, condc, bnrc, CHI;
-   if ((flag < 0) || (flag > 1))
+   if (flag < 0 || flag > 1) err(flagerr,"bnrrootnumber");
-     err(flagerr,"bnrrootnumber");
    checkbnr(bnr);
    cond = gmael(bnr, 2, 1);
    l    = lg(gmael(bnr, 5, 2));
-   Pi   = mppi(prec);
-   if ((typ(chi) != t_VEC) || (lg(chi) != l))
+   if (typ(chi) != t_VEC || lg(chi) != l)
      err(talker, "incorrect character in bnrrootnumber");
-   if (!flag)
+   if (flag) condc = cond;
+   else
    {
      condc = bnrconductorofchar(bnr, chi);
      if (gegal(cond, condc)) flag = 1;
    }
-   else condc = cond;
    if (flag)
+   {
+     GEN cyc  = gmael(bnr, 5, 2);
      bnrc = bnr;
+     CHI = get_Char(get_chic(chi,cyc), prec);
+   }
    else
+   {
      bnrc = buchrayinitgen((GEN)bnr[1], condc);
+     CHI = (GEN)GetPrimChar(chi, bnr, bnrc, prec)[1];
-   chic = cgetg(l, t_VEC);
+   }
-   cyc  = gmael(bnr, 5, 2);
+   return gerepilecopy(av, (GEN)ComputeArtinNumber(bnrc, _vec(CHI), 1, prec)[1]);
-   for (i = 1; i < l; i++)
-     chic[i] = ldiv((GEN)chi[i], (GEN)cyc[i]);
-   d = denom(chic);
-   p2 = cgetg(4, t_VEC);
-   p2[1] = lmul(d, chic);
-   if (egalii(d, gdeux))
-     p2[2] = lstoi(-1);
-   else
-     p2[2] = lexp(gdiv(gmul2n(gmul(gi, Pi), 1), d), prec);
-   p2[3] = (long)d;
-   dtcr = cgetg(9, t_VEC);
-   dtcr[1] = (long)chi;
-   dtcr[2] = zero;
-   dtcr[3] = (long)bnrc;
-   dtcr[4] = (long)bnr;
-   dtcr[5] = (long)p2;
-   dtcr[6] = zero;
-   dtcr[7] = (long)condc;
-   p1 = GetPrimChar(chi, bnr, bnrc, prec);
-   if (!p1)
-     dtcr[8] = (long)p2;
-   else
-     dtcr[8] = p1[1];
-   return gerepileupto(av, ComputeArtinNumber(dtcr, 0, prec));
  }
  /********************************************************************/
-Line 909  bnrrootnumber(GEN bnr, GEN chi, long flag, long prec)
+Line 919  bnrrootnumber(GEN bnr, GEN chi, long flag, long prec)
 Line 909  bnrrootnumber(GEN bnr, GEN chi, long flag, long prec)
 Line 919  bnrrootnumber(GEN bnr, GEN chi, long flag, long prec)
  static GEN
  LiftChar(GEN cyc, GEN Mat, GEN chi)
  {
-   long lm, l, i, j, av;
+   long lm, l, i, j;
+   gpmem_t av;
    GEN lchi, s;
    lm = lg(cyc) - 1;
-Line 936  LiftChar(GEN cyc, GEN Mat, GEN chi)
+Line 947  LiftChar(GEN cyc, GEN Mat, GEN chi)
 Line 936  LiftChar(GEN cyc, GEN Mat, GEN chi)
 Line 947  LiftChar(GEN cyc, GEN Mat, GEN chi)
  static GEN
  ComputeAChi(GEN dtcr, long flag, long prec)
  {
-   long l, i;
+   long l, i, r;
-   GEN p1, ray, r, A, rep, diff, chi, bnrc;
+   GEN p1, ray, A, rep, diff, chi, bnrc, nf;
    diff = (GEN)dtcr[6];
-   bnrc = (GEN)dtcr[3];
+   bnrc = (GEN)dtcr[3]; nf = checknf(bnrc);
    chi  = (GEN)dtcr[8];
    l    = lg(diff) - 1;
    A = gun;
-   r = gzero;
+   r = 0;
    for (i = 1; i <= l; i++)
    {
+     GEN B;
      ray = isprincipalray(bnrc, (GEN)diff[i]);
-     p1  = ComputeImagebyChar(chi, ray, 0);
+     p1  = ComputeImagebyChar(chi, ray);
      if (flag)
-       A = gmul(A, gsub(gun, gdiv(p1, idealnorm((GEN)bnrc[1], (GEN)diff[i]))));
+       B = gsub(gun, gdiv(p1, idealnorm(nf, (GEN)diff[i])));
-     else
+     else if (gcmp1(p1))
      {
-       if (gcmp1(p1))
+       B = glog(idealnorm(nf, (GEN)diff[i]), prec);
-       {
+       r++;
-         r = addis(r, 1);
-         A = gmul(A, glog(idealnorm((GEN)bnrc[1], (GEN)diff[i]), prec));
-       }
-       else
-         A = gmul(A, gsub(gun, p1));
      }
+     else
+       B = gsub(gun, p1);
+     A = gmul(A, B);
    }
    if (flag) return A;
    rep = cgetg(3, t_VEC);
-   rep[1] = (long)r;
+   rep[1] = lstoi(r);
    rep[2] = (long)A;
    return rep;
  }
+ static GEN
+ _data9(GEN arch, long r1, long r2)
+ {
+   GEN z = cgetg(5, t_VECSMALL);
+   long i, b, q = 0;
+   for (i=1; i<=r1; i++) if (signe(arch[i])) q++;
+   z[1] = q; b = r1 - q;
+   z[2] = b;
+   z[3] = r2;
+   z[4] = max(b+r2+1, r2+q);
+   return z;
+ }
  /* Given a list [chi, cond(chi)] of characters over Cl(bnr), compute a
     vector dataCR containing for each character:
 : chi
-Line 992  static GEN
+Line 1015  static GEN
 Line 992  static GEN
 Line 1015  static GEN
  InitChar(GEN bnr, GEN listCR, long prec)
  {
    GEN bnf = checkbnf(bnr), nf = checknf(bnf);
-   GEN modul, dk, C, dataCR, chi, cond, Mr, chic, gr2, p1, p2, q;
+   GEN modul, dk, C, dataCR, chi, cond, Mr, z1, p1;
-   long N, r1, r2, prec2, h, i, j, nbg, av = avma;
+   long N, r1, r2, prec2, h, i, j;
+   gpmem_t av = avma;
    modul = gmael(bnr, 2, 1);
    Mr    = gmael(bnr, 5, 2);
-   nbg   = lg(Mr) - 1;
    dk    = (GEN)nf[3];
    N     = degpol(nf[1]);
-   r1    = nf_get_r1(nf);
+   nf_get_sign(nf, &r1,&r2);
-   r2    = (N - r1) >> 1; gr2 = stoi(r2);
+   prec2 = ((prec-2) << 1) + EXTRA_PREC;
-   prec2 = ((prec - 2)<<1) + EXTRA_PREC;
    C     = gmul2n(gsqrt(gdiv(absi(dk), gpowgs(mppi(prec2),N)), prec2), -r2);
    disable_dbg(0);
-Line 1012  InitChar(GEN bnr, GEN listCR, long prec)
+Line 1034  InitChar(GEN bnr, GEN listCR, long prec)
 Line 1012  InitChar(GEN bnr, GEN listCR, long prec)
 Line 1034  InitChar(GEN bnr, GEN listCR, long prec)
    for (i = 1; i <= h ;i++)
      dataCR[i] = lgetg(10, t_VEC);
-   q = gnorml2((GEN)modul[2]);
+   z1 = _data9((GEN)modul[2],r1,r2);
-   p1 = cgetg(5, t_VEC);
-   p1[1] = (long)q;
-   p1[2] = lsubsi(r1, q);
-   p1[3] = (long)gr2;
-   p1[4] = lmax(addis((GEN)p1[2], r2+1), addsi(r2, q));
    for (i = 1; i <= h; i++)
    {
-Line 1039  InitChar(GEN bnr, GEN listCR, long prec)
+Line 1056  InitChar(GEN bnr, GEN listCR, long prec)
 Line 1039  InitChar(GEN bnr, GEN listCR, long prec)
 Line 1056  InitChar(GEN bnr, GEN listCR, long prec)
        data[3] = (long)bnr;
        data[6] = lgetg(1, t_VEC);
        data[7] = modul[1];
-       data[9] = (long)p1;
+       data[9] = (long)z1;
        olddata = data;
      }
-Line 1056  InitChar(GEN bnr, GEN listCR, long prec)
+Line 1073  InitChar(GEN bnr, GEN listCR, long prec)
 Line 1056  InitChar(GEN bnr, GEN listCR, long prec)
 Line 1073  InitChar(GEN bnr, GEN listCR, long prec)
        data[3] = olddata[3]; /* bnr(cond(chi)) */
      }
      data[4] = (long)bnr; /* bnr(m) */
+     data[5] = (long)get_Char(get_chic(chi,Mr),prec); /* associated to bnr(m) */
-     chic = cgetg(nbg + 1, t_VEC);
-     for (j = 1; j <= nbg; j++)
-       chic[j] = ldiv((GEN)chi[j], (GEN)Mr[j]);
-     data[5] = (long)get_Char(chic,prec); /* char associated to bnr(m) */
      /* compute diff(chi) and the corresponding primitive character */
      data[7] = cond[1];
-     p2 = GetPrimChar(chi, bnr, (GEN)data[3], prec2);
+     p1 = GetPrimChar(chi, bnr, (GEN)data[3], prec2);
-     if (p2)
+     if (p1)
      {
-       data[6] = p2[2];
+       data[6] = p1[2];
-       data[8] = p2[1];
+       data[8] = p1[1];
      }
      else
      {
        data[6] = lgetg(1, t_VEC);
        data[8] = data[5];
      }
+     data[9] = olddata? olddata[9]: (long)_data9((GEN)cond[2],r1,r2);
-     /* compute q and store [q, r1 - q, r2] */
-     if (!olddata)
-     {
-       q = gnorml2((GEN)cond[2]);
-       p2 = cgetg(5, t_VEC);
-       p2[1] = (long)q;
-       p2[2] = lsubsi(r1, q);
-       p2[3] = (long)gr2;
-       p2[4] = lmax(addis((GEN)p2[2], r2+1), addsi(r2, q));
-       data[9] = (long)p2;
-     }
-     else
-       data[9] = olddata[9];
    }
    disable_dbg(-1);
 Line 1101  static GEN
 Line 1101  static GEN
 Line 1101  static GEN
  InitChar0(GEN dataD, long prec)
  {
    GEN MrD, listCR, p1, chi, lchi, Surj, cond, bnr, p2, Mr, d, allCR;
-   long hD, h, nc, i, j, lD, nbg, tnc, av = avma;
+   long hD, h, nc, i, j, lD, tnc;
+   gpmem_t av = avma;
    Surj = gmael(dataD, 2, 3);
    MrD  = gmael(dataD, 2, 2);
-Line 1110  InitChar0(GEN dataD, long prec)
+Line 1111  InitChar0(GEN dataD, long prec)
 Line 1110  InitChar0(GEN dataD, long prec)
 Line 1111  InitChar0(GEN dataD, long prec)
    hD   = itos(gmael(dataD, 2, 1));
    h    = hD >> 1;
    lD   = lg(MrD)-1;
-   nbg  = lg(Mr) - 1;
    disable_dbg(0);
 Line 1119  InitChar0(GEN dataD, long prec)
 Line 1119  InitChar0(GEN dataD, long prec)
 Line 1119  InitChar0(GEN dataD, long prec)
    allCR  = cgetg(h + 1, t_VEC); /* all characters, including conjugates */
    tnc = 1;
-   p1 = FindEltofGroup(hD, MrD);
+   p1 = EltsOfGroup(hD, MrD);
    for (i = 1; tnc <= h; i++)
    {
 Line 1142  InitChar0(GEN dataD, long prec)
 Line 1142  InitChar0(GEN dataD, long prec)
 Line 1142  InitChar0(GEN dataD, long prec)
      listCR[nc++] = (long)p2;
      allCR[tnc++] = (long)lchi;
-     /* compute the order of chi */
-     p2 = cgetg(nbg + 1, t_VEC);
-     for (j = 1; j <= nbg; j++)
-       p2[j] = ldiv((GEN)lchi[j], (GEN)Mr[j]);
-     d = denom(p2);
      /* if chi is not real, add its conjugate character to allCR */
-     if (!gegal(d, gdeux))
+     d = Order(Mr, lchi);
+     if (!egalii(d, gdeux))
        allCR[tnc++] = (long)ConjChar(lchi, Mr);
    }
-Line 1162  InitChar0(GEN dataD, long prec)
+Line 1157  InitChar0(GEN dataD, long prec)
 Line 1162  InitChar0(GEN dataD, long prec)
 Line 1157  InitChar0(GEN dataD, long prec)
  static GEN
  CharNewPrec(GEN dataCR, GEN nf, long prec)
  {
-   GEN dk, C, p1, Pi;
+   GEN dk, C, p1;
-   long av = avma, N, l, j, prec2;
+   long N, l, j, prec2;
    dk    =  (GEN)nf[3];
    N     =  degpol((GEN)nf[1]);
-   l     =  lg(dataCR) - 1;
    prec2 = ((prec - 2)<<1) + EXTRA_PREC;
-   Pi    = mppi(prec2);
+   C = mpsqrt(gdiv(dk, gpowgs(mppi(prec2), N)));
-   C = gsqrt(gdiv(dk, gpowgs(Pi, N)), prec2);
+   l = lg(dataCR);
+   for (j = 1; j < l; j++)
-   for (j = 1; j <= l; j++)
    {
-     mael(dataCR, j, 2) = lmul(C, gsqrt(det(gmael(dataCR, j, 7)), prec2));
+     GEN dtcr = (GEN)dataCR[j];
+     dtcr[2] = lmul(C, gsqrt(dethnf_i((GEN)dtcr[7]), prec2));
-     mael4(dataCR, j, 3, 1, 7) = lcopy(nf);
+     mael3(dtcr, 3, 1, 7) = (long)nf;
-     p1 = gmael(dataCR, j, 5);
+     p1 = (GEN)dtcr[5]; p1[2] = (long)InitRU((GEN)p1[3], prec2);
-     p1[2] = lexp(gdiv(gmul2n(gmul(gi, Pi), 1), (GEN)p1[3]), prec);
+     p1 = (GEN)dtcr[8]; p1[2] = (long)InitRU((GEN)p1[3], prec2);
-     p1 = gmael(dataCR, j, 8);
-     p1[2] = lexp(gdiv(gmul2n(gmul(gi, Pi), 1), (GEN)p1[3]), prec);
    }
-   return gerepilecopy(av, dataCR);
+   return dataCR;
  }
  /********************************************************************/
-Line 1198  CharNewPrec(GEN dataCR, GEN nf, long prec)
+Line 1189  CharNewPrec(GEN dataCR, GEN nf, long prec)
 Line 1198  CharNewPrec(GEN dataCR, GEN nf, long prec)
 Line 1189  CharNewPrec(GEN dataCR, GEN nf, long prec)
  /********************************************************************/
  static void
- _0toCoeff(int *rep, long dg)
+ _0toCoeff(int *rep, long deg)
  {
    long i;
-   for (i=0; i<dg; i++) rep[i] = 0;
+   for (i=0; i<deg; i++) rep[i] = 0;
  }
  /* transform a polmod into coeff */
  static void
- Polmod2Coeff(int *rep, GEN polmod, long dg)
+ Polmod2Coeff(int *rep, GEN polmod, long deg)
  {
    GEN pol = (GEN)polmod[2];
    long i,d = degpol(pol);
    pol += 2;
    for (i=0; i<=d; i++) rep[i] = itos((GEN)pol[i]);
-   for (   ; i<dg; i++) rep[i] = 0;
+   for (   ; i<deg; i++) rep[i] = 0;
  }
- /* initialize a cl x nmax x [degs[1], ..., degs[cl] matrix of ints */
+ /* initialize a deg * n matrix of ints */
- /* modified to allocate ints and pointers separately since this used to
+ static int**
-    break on 64bit platforms --GN1999Sep01 */
+ InitMatAn(long n, long deg, long flag)
- static int***
- InitMatAn(long cl, long nmax, GEN degs, long flag)
  {
-   long si,dg,i,j,k, n = nmax+1;
+   long i, j;
-   int *c, **pj, ***A;
+   int *a, **A = (int**)gpmalloc((n+1)*sizeof(int*));
+   A[0] = NULL;
-   for (si=0, i=1; i <= cl; i++)
+   for (i = 1; i <= n; i++)
    {
-     dg = degs[i];
+     a = (int*)gpmalloc(deg*sizeof(int));
-     si += dg;
+     A[i] = a; a[0] = (i == 1 || flag);
+     for (j = 1; j < deg; j++) a[j] = 0;
    }
-   A = (int***)gpmalloc((cl+1)*sizeof(int**) + cl*n*sizeof(int*));
-   c = (int*)gpmalloc(si*n*sizeof(int));
-   A[0] = (int**)c;              /* keep it around for FreeMat() */
-   pj = (int**)(A + (cl+1));
-   for (i = 1; i <= cl; i++, pj+=n)
-   {
-     A[i] = pj; dg = degs[i];
-     for (j = 0; j < n; j++,c+=dg)
-     {
-       pj[j] = c;
-       c[0] = (j == 1 || flag);
-       for (k = 1; k < dg; k++) c[k] = 0;
-     }
-   }
    return A;
  }
  static void
- FreeMat(int ***m)
+ FreeMat(int **A, long n)
  {
-   free((void *)(m[0]));
+   long i;
-   free((void *)m);
+   for (i = 0; i <= n; i++)
+     if (A[i]) free((void*)A[i]);
+   free((void*)A);
  }
  /* initialize coeff reduction */
- /* similar change --GN1999Sep01 */
+ static int**
- static int***
+ InitReduction(GEN CHI, long deg)
- InitReduction(GEN dataCR, GEN degs)
  {
-   long av = avma,si,sp,dg,i,j, cl = lg(dataCR)-1;
+   long j;
-   int *c, **pj, ***A;
+   gpmem_t av = avma;
+   int **A;
    GEN d,polmod,pol, x = polx[0];
-   for (si=sp=0, i=1; i <= cl; i++)
+   A   = (int**)gpmalloc(deg*sizeof(int*));
+   d   = (GEN)CHI[3];
+   pol = cyclo(itos(d), 0);
+   for (j = 0; j < deg; j++)
    {
-     dg = degs[i];
+     A[j] = (int*)gpmalloc(deg*sizeof(int));
-     sp += dg;
+     polmod = gmodulcp(gpowgs(x, deg + j), pol);
-     si += dg*dg;
+     Polmod2Coeff(A[j], polmod, deg);
    }
-   A = (int***)gpmalloc((cl+1)*sizeof(int**) + sp*sizeof(int*));
-   c = (int*)gpmalloc(si*sizeof(int));
-   A[0] = (int**)c;              /* keep it around for FreeMat() */
-   pj = (int**)(A + (cl+1));
-   for (i = 1; i <= cl; i++)
-   {
-     A[i] = pj;
-     d   = gmael3(dataCR, i, 5, 3);
-     pol = cyclo(itos(d), 0);
-     dg  = degs[i]; /* degree(pol) */
-     for (j = 0; j < dg; j++,c+=dg)
-     {
-       pj[j] = c;
-       polmod = gmodulcp(gpowgs(x, dg + j), pol);
-       Polmod2Coeff(c, polmod, dg);
-     }
-     pj += dg;
-   }
    avma = av; return A;
  }
  #if 0
- static void
+ void
- pan(int ***an,long cl, long nmax, GEN degs)
+ pan(int **an, long n, long deg)
  {
-   long i,j,k;
+   long i,j;
-   for (i = 1; i <= cl; i++)
+   for (i = 1; i <= n; i++)
    {
-     long dg = degs[i];
+     fprintferr("n = %ld: ",i);
-     for (j = 0; j <= nmax; j++)
+     for (j = 0; j < deg; j++) fprintferr("%d ",an[i][j]);
-       for (k = 0; k < dg; k++) fprintferr("%d ",an[i][j][k]);
+     fprintferr("\n");
    }
  }
  #endif
- /* multiply (with reduction) a polmod with a coeff. put result in c1 */
+ /* returns 0 if c is zero, 1 otherwise. */
+ static int
+ IsZero(int* c, long deg)
+ {
+   long i;
+   for (i = 0; i < deg; i++)
+     if (c[i]) return 0;
+   return 1;
+ }
+ /* set c0 <-- c0 * c1 */
  static void
- MulPolmodCoeff(GEN polmod, int* c1, int** reduc, long dg)
+ MulCoeff(int *c0, int* c1, int** reduc, long deg)
  {
-   long av,i,j;
+   long i,j;
-   int c, *c2, *c3;
+   int c, *T;
-   if (gcmp1(polmod)) return;
+   if (IsZero(c0,deg)) return;
-   for (i = 0; i < dg; i++)
-     if (c1[i]) break;
-   if (i == dg) return;
-   av = avma;
-   c3 = (int*)new_chunk(2*dg);
-   c2 = (int*)new_chunk(dg);
-   Polmod2Coeff(c2,polmod, dg);
-   for (i = 0; i < 2*dg; i++)
+   T = (int*)new_chunk(2*deg);
+   for (i = 0; i < 2*deg; i++)
    {
      c = 0;
      for (j = 0; j <= i; j++)
-       if (j < dg && j > i - dg) c += c1[j] * c2[i-j];
+       if (j < deg && j > i - deg) c += c0[j] * c1[i-j];
-     c3[i] = c;
+     T[i] = c;
    }
-   c2 = c1;
+   for (i = 0; i < deg; i++)
-   for (i = 0; i < dg; i++)
    {
-     c = c3[i];
+     c = T[i];
-     for (j = 0; j < dg; j++) c += reduc[j][i] * c3[dg+j];
+     for (j = 0; j < deg; j++) c += reduc[j][i] * T[deg+j];
-     c2[i] = c;
+     c0[i] = c;
    }
-   /* cast necessary to work around a gcc-2.96 bug on alpha-linux (IS) */
-   for (     ; i < (short)dg; i++) c2[i] = 0;
-   avma = av;
  }
- /* c0 <- c0 + c2 * c1 */
+ /* c0 <- c0 + c1 * c2 */
  static void
- AddMulCoeff(int *c0, int *c2, int* c1, int** reduc, long dg)
+ AddMulCoeff(int *c0, int *c1, int* c2, int** reduc, long deg)
  {
-   long av,i,j;
+   long i, j;
-   int c, *c3;
+   gpmem_t av;
+   int c, *t;
-   if (!c2) /* c2 == 1 */
+   if (IsZero(c2,deg)) return;
+   if (!c1) /* c1 == 1 */
    {
-     for (i = 0; i < dg; i++) c0[i] += c1[i];
+     for (i = 0; i < deg; i++) c0[i] += c2[i];
      return;
    }
-   for (i = 0; i <= dg; i++)
-     if (c1[i]) break;
-   if (i > dg) return;
    av = avma;
-   c3 = (int*)new_chunk(2*dg);
+   t = (int*)new_chunk(2*deg); /* = c1 * c2, not reduced */
-   for (i = 0; i < 2*dg; i++)
+   for (i = 0; i < 2*deg; i++)
    {
      c = 0;
      for (j = 0; j <= i; j++)
-       if (j < dg && j > i - dg) c += c1[j] * c2[i-j];
+       if (j < deg && j > i - deg) c += c1[j] * c2[i-j];
-     c3[i] = c;
+     t[i] = c;
    }
-   for (i = 0; i < dg; i++)
+   for (i = 0; i < deg; i++)
    {
-     c = c0[i] + c3[i];
+     c = t[i];
-     for (j = 0; j < dg; j++) c += reduc[j][i] * c3[dg+j];
+     for (j = 0; j < deg; j++) c += reduc[j][i] * t[deg+j];
-     c0[i] = c;
+     c0[i] += c;
    }
    avma = av;
  }
- /* returns 0 if c is zero, 1 otherwise. */
- static long
- IsZero(int* c, long dg)
- {
-   long i;
-   for (i = 0; i < dg; i++)
-     if (c[i]) return 0;
-   return 1;
- }
  /* evaluate the coeff. No Garbage collector */
  static GEN
- EvalCoeff(GEN z, int* c, long dg)
+ EvalCoeff(GEN z, int* c, long deg)
  {
    long i,j;
    GEN e, r;
+   if (!c) return gzero;
  #if 0
    /* standard Horner */
-   e = stoi(c[dg - 1]);
+   e = stoi(c[deg - 1]);
-   for (i = dg - 2; i >= 0; i--)
+   for (i = deg - 2; i >= 0; i--)
      e = gadd(stoi(c[i]), gmul(z, e));
  #else
    /* specific attention to sparse polynomials */
    e = NULL;
-   for (i = dg-1; i >=0; i=j-1)
+   for (i = deg-1; i >=0; i=j-1)
    {
      for (j=i; c[j] == 0; j--)
        if (j==0)
-Line 1422  EvalCoeff(GEN z, int* c, long dg)
+Line 1372  EvalCoeff(GEN z, int* c, long dg)
 Line 1422  EvalCoeff(GEN z, int* c, long dg)
 Line 1372  EvalCoeff(GEN z, int* c, long dg)
    return e;
  }
- /* copy the n * m array matan */
+ /* copy the n * (m+1) array matan */
  static void
- CopyCoeff(int*** a, int*** a2, long n, long m, GEN degs)
+ CopyCoeff(int** a, int** a2, long n, long m)
  {
-   long i,j,k;
+   long i,j;
    for (i = 1; i <= n; i++)
    {
-     long dg = degs[i];
+     int *b = a[i], *b2 = a2[i];
-     int **b = a[i], **b2 = a2[i];
+     for (j = 0; j < m; j++) b2[j] = b[j];
-     for (j = 0; j <= m; j++)
-     {
-       int *c = b[j], *c2 = b2[j];
-       for (k = 0; k < dg; k++) c2[k] = c[k];
-     }
    }
-   return;
  }
- /* initialize the data for GetRay */
+ /* return q*p if <= n. Beware overflow */
- static GEN
+ static long
- InitGetRay(GEN bnr,  long nmax)
+ next_pow(long q, long p, long n)
  {
-   long bd, i, j, l;
+   const GEN x = muluu((ulong)q, (ulong)p);
-   GEN listid, listcl, id, rep, bnf, cond;
+   const ulong qp = (ulong)x[2];
+   return (lgefint(x) > 3 || qp > (ulong)n)? 0: qp;
-   bnf  =  (GEN)bnr[1];
-   cond =  gmael3(bnr, 2, 1, 1);
-   if (nmax < 1000) return NULL;
-   rep = cgetg(4, t_VEC);
-   disable_dbg(0);
-   bd = min(1000, nmax / 50);
-   listid = ideallist(bnf, bd);
-   disable_dbg(-1);
-   listcl = cgetg(bd + 1, t_VEC);
-   for (i = 1; i <= bd; i++)
-   {
-     l = lg((GEN)listid[i]) - 1;
-     listcl[i] = lgetg(l + 1, t_VEC);
-     for (j = 1; j <= l; j++)
-     {
-       id = gmael(listid, i, j);
-       if (gcmp1(gcoeff(idealadd(bnf, id, cond), 1, 1)))
-         mael(listcl, i, j) = (long)isprincipalray(bnr, id);
-     }
-   }
-   if (DEBUGLEVEL) msgtimer("InitGetRay");
-   rep[1] = (long)listid;
-   rep[2] = (long)listcl;
-   rep[3] = nf_get_r2((GEN)bnf[7])? 0: un; /* != 0 iff nf is totally real */
-   return rep;
  }
- /* compute the class of the prime ideal pr in cl(bnr) using dataray */
+ static void
- static GEN
+ an_AddMul(int **an,int **an2, long np, long n, long deg, GEN chi, int **reduc)
- GetRay(GEN bnr,  GEN dataray,  GEN pr, long prec)
  {
-   long av = avma, N, n, bd, c;
+   GEN chi2 = chi;
-   GEN id, tid, t2, u, alpha, p1, cl, listid, listcl, nf, cond;
+   long q, qk, k;
+   int *c, *c2 = (int*)new_chunk(deg);
-   if (!dataray)
+   CopyCoeff(an, an2, n/np, deg);
-     return isprincipalray(bnr, pr);
+   for (q=np;;)
-   listid =  (GEN)dataray[1];
-   listcl =  (GEN)dataray[2];
-   cond   =  gmael3(bnr, 2, 1, 1);
-   bd     =  lg(listid) - 1;
-   nf     =  gmael(bnr, 1, 7);
-   N      =  degpol(nf[1]);
-   if (dataray[3])
-     t2 = gmael(nf, 5, 4);
-   else
-     t2 = gmael(nf, 5, 3);
-   id  = prime_to_ideal(nf, pr);
-   tid = qf_base_change(t2, id, 1);
-   if (dataray[3])
-     u = lllgramint(tid);
-   else
-     u = lllgramintern(tid,100,1,prec);
-   if (!u) return gerepileupto(av, isprincipalray(bnr, id));
-   c = 1; alpha = NULL;
-   for (c=1; c<=N; c++)
    {
-     p1 = gmul(id, (GEN)u[c]);
+     if (gcmp1(chi2)) c = NULL; else { Polmod2Coeff(c2, chi2, deg); c = c2; }
-     if (gcmp1(gcoeff(idealadd(nf, p1, cond), 1, 1))) { alpha = p1; break; }
+     for(k = 1, qk = q; qk <= n; k++, qk += q)
-   }
+       AddMulCoeff(an[qk], c, an2[k], reduc, deg);
-   if (!alpha)
+     if (! (q = next_pow(q,np, n)) ) break;
-     return gerepileupto(av, isprincipalray(bnr, pr));
-   id = idealdivexact(nf, alpha, id);
+     chi2 = gmul(chi2, chi);
-   n = itos(det(id));
-   if (n > bd)
-     cl = isprincipalray(bnr, id);
-   else
-   {
-     cl = NULL;
-     c  = 1;
-     p1 = (GEN)listid[n];
-     while (!cl)
-     {
-       if (gegal((GEN)p1[c], id))
-         cl = gmael(listcl, n, c);
-       c++;
-     }
    }
-   return gerepileupto(av, gsub(isprincipalray(bnr, alpha), cl));
  }
  /* correct the coefficients an(chi) according with diff(chi) in place */
  static void
- CorrectCoeff(GEN dtcr, int** an, int** reduc, long nmax, long dg)
+ CorrectCoeff(GEN dtcr, int** an, int** reduc, long n, long deg)
  {
-   long lg, av1, j, p, q, limk, k, l, av = avma;
+   gpmem_t av = avma;
-   int ***an2, **an1, *c, *c2;
+   long lg, j, np;
-   GEN chi, bnrc, diff, ray, ki, ki2, pr, degs;
+   gpmem_t av1;
+   int **an2;
+   GEN bnrc, diff, ray, chi, CHI, pr;
+   CHI_t C;
-   chi  =  (GEN)dtcr[8];
-   bnrc =  (GEN)dtcr[3];
    diff =  (GEN)dtcr[6];
    lg   =  lg(diff) - 1;
    if (!lg) return;
-   if (DEBUGLEVEL > 2) fprintferr("diff(chi) = %Z", diff);
+   bnrc =  (GEN)dtcr[3];
+   CHI  =  (GEN)dtcr[8]; init_CHI_alg(&C, CHI);
-   degs = cgetg(2, t_VECSMALL); degs[1] = dg;
+   if (DEBUGLEVEL > 2) fprintferr("diff(CHI) = %Z", diff);
-   an2 = InitMatAn(1, nmax, degs, 0); an1 = an2[1];
-   c = (int*)new_chunk(dg);
-   av1 = avma;
+   an2 = InitMatAn(n, deg, 0);
+   av1 = avma;
    for (j = 1; j <= lg; j++)
    {
-     for (k = 0; k <= nmax; k++)
-       for (l = 0; l < dg; l++) an1[k][l] = an[k][l];
      pr  = (GEN)diff[j];
+     np = itos(powgi((GEN)pr[1], (GEN)pr[4]));
      ray = isprincipalray(bnrc, pr);
-     ki  = ComputeImagebyChar(chi, ray, 1);
+     chi  = EvalChar(&C,ray);
-     ki2 = gcopy(ki);
-     q = p = itos(powgi((GEN)pr[1], (GEN)pr[4]));
+     an_AddMul(an,an2,np,n,deg,chi,reduc);
-     limk = nmax / q;
-     while (q <= nmax)
-     {
-       if (gcmp1(ki2)) c2 = NULL; else { Polmod2Coeff(c,ki2, dg); c2 = c; }
-       for(k = 1; k <= limk; k++)
-         AddMulCoeff(an[k*q], c2, an1[k], reduc, dg);
-       q *= p; limk /= p;
-       ki2 = gmul(ki2, ki);
-     }
      avma = av1;
    }
-   FreeMat(an2); avma = av;
+   FreeMat(an2, n); avma = av;
  }
  /* compute the coefficients an in the general case */
- static int***
+ static int**
- ComputeCoeff(GEN dataCR, long nmax, long prec)
+ ComputeCoeff(GEN dtcr, LISTray *R, long n, long deg)
  {
-   long cl, i, j, av = avma, av2, np, q, q1, limk, k, id, cpt = 10, dg, Bq;
+   gpmem_t av = avma, av2;
-   int ***matan, ***reduc, ***matan2, *c2;
+   long i, l, np;
-   GEN c, degs, tabprem, bnf, pr, cond, ray, ki, ki2, prime, npg, bnr, dataray;
+   int **an, **reduc, **an2;
-   byteptr dp = diffptr;
+   GEN L, CHI, ray, chi;
+   CHI_t C;
-   bnr  =  gmael(dataCR, 1, 4);
+   CHI = (GEN)dtcr[5]; init_CHI_alg(&C, CHI);
-   bnf  =  (GEN)bnr[1];
-   cond =  gmael3(bnr, 2, 1, 1);
-   cl   =  lg(dataCR) - 1;
-   dataray = InitGetRay(bnr, nmax);
+   an  = InitMatAn(n, deg, 0);
+   an2 = InitMatAn(n, deg, 0);
+   reduc  = InitReduction(CHI, deg);
+   av2 = avma;
-   degs = GetDeg(dataCR);
+   L = R->L1; l = lg(L);
-   matan  = InitMatAn(cl, nmax, degs, 0);
+   for (i=1; i<l; i++, avma = av2)
-   matan2 = InitMatAn(cl, nmax, degs, 0);
+   {
-   reduc  = InitReduction(dataCR, degs);
+     np = L[i]; ray = R->L1ray[i];
-   c = cgetg(cl + 1, t_VEC);
+     chi  = EvalChar(&C, ray);
-   for (i = 1; i <= cl; i++)
+     an_AddMul(an,an2,np,n,deg,chi,reduc);
-     c[i] = (long)new_chunk(degs[i]);
+   }
+   FreeMat(an2, n);
-   if (DEBUGLEVEL > 1) fprintferr("p = ");
+   CorrectCoeff(dtcr, an, reduc, n, deg);
+   FreeMat(reduc, deg-1);
+   avma = av; return an;
+ }
-   prime = stoi(2); dp++;
+ /********************************************************************/
-   av2 = avma;
+ /*              5th part: compute L-functions at s=1                */
-   while (*dp && (prime[2] <= nmax))
+ /********************************************************************/
-   {
+ static void
-     tabprem = primedec(bnf, prime);
+ deg11(LISTray *R, long p, GEN bnr, GEN pr) {
-     for (j = 1; j < lg(tabprem); j++)
+   GEN z = isprincipalray(bnr, pr);
-     {
+   appendL(R->L1, (GEN)p);
-       pr  = (GEN)tabprem[j];
+   appendL((GEN)R->L1ray, z);
-       npg = powgi((GEN)pr[1], (GEN)pr[4]);
+ }
-       if (is_bigint(npg) || (np=npg[2]) > nmax
+ static void
-                          || idealval(bnf, cond, pr)) continue;
+ deg12(LISTray *R, long p, GEN bnr, GEN Lpr) {
+   GEN z = isprincipalray(bnr, (GEN)Lpr[1]);
+   appendL(R->L11, (GEN)p);
+   appendL((GEN)R->L11ray, z);
+ }
+ static void
+ deg0(LISTray *R, long p) {
+   appendL(R->L0, (GEN)p);
+ }
+ static void
+ deg2(LISTray *R, long p) {
+   appendL(R->L2, (GEN)p);
+ }
-       CopyCoeff(matan, matan2, cl, nmax, degs);
+ /* pi(x) <= ?? */
-       ray = GetRay(bnr, dataray, pr, prec);
+ static long
-       ki  = chiideal(dataCR, ray, 1);
+ PiBound(long x)
-       ki2 = dummycopy(ki);
+ {
+   double lx = log((double)x);
+   return 1 + (long) (x/lx * (1 + 3/(2*lx)));
+ }
-       Bq = nmax/np;
+ static void
-       for (q1 = 1; q1 <= Bq; q1 *= np)
+ InitPrimesQuad(GEN bnr, long nmax, LISTray *R)
+ {
+   gpmem_t av = avma;
+   GEN bnf = (GEN)bnr[1], cond = gmael3(bnr,2,1,1);
+   long p,i,l, condZ = itos(gcoeff(cond,1,1)), contZ = itos(content(cond));
+   GEN prime, pr, nf = checknf(bnf), dk = (GEN)nf[3];
+   byteptr d = diffptr + 1;
+   GEN *gptr[7];
+   l = 1 + PiBound(nmax);
+   R->L0  = cget1(l, t_VECSMALL);
+   R->L2  = cget1(l, t_VECSMALL); R->condZ = condZ;
+   R->L1 = cget1(l, t_VECSMALL); R->L1ray = (GEN*)cget1(l, t_VEC);
+   R->L11 = cget1(l, t_VECSMALL); R->L11ray = (GEN*)cget1(l, t_VEC);
+   prime = stoi(2);
+   for (p = 2; p <= nmax; prime[2] = p) {
+     switch (krogs(dk, p))
+     {
+     case -1: /* inert */
+       if (condZ % p == 0) deg0(R,p); else deg2(R,p);
+       break;
+     case 1: /* split */
+       pr = primedec(nf, prime);
+       if      (condZ % p != 0) deg12(R, p, bnr, pr);
+       else if (contZ % p == 0) deg0(R,p);
+       else
        {
-         q = q1*np;
+         pr = idealval(nf, cond, (GEN)pr[1])? (GEN)pr[2]: (GEN)pr[1];
-         limk = nmax / q;
+         deg11(R, p, bnr, pr);
-         for (id = 1; id <= cl; id++)
-         {
-           dg = degs[id];
-           if (gcmp1((GEN)ki2[id]))
-             c2 = NULL;
-           else
-           {
-             c2 = (int*)c[id];
-             Polmod2Coeff(c2, (GEN)ki2[id], dg);
-           }
-           for (k = 1; k <= limk; k++)
-             AddMulCoeff(matan[id][k*q], c2, matan2[id][k], reduc[id], dg);
-           ki2[id] = lmul((GEN)ki2[id], (GEN)ki[id]);
-         }
        }
+       break;
+     default: /* ramified */
+       if (condZ % p == 0) deg0(R,p);
+       else
+       {
+         pr = (GEN)primedec(nf, prime)[1];
+         deg11(R, p, bnr, pr);
+       }
+       break;
      }
-     avma = av2;
+     NEXT_PRIME_VIADIFF(p,d);
-     prime[2] += (*dp++);
+   }
-     if (!*dp) err(primer1);
+   /* precompute isprincipalray(x), x in Z */
+   R->rayZ = (GEN*)cgetg(condZ, t_VEC);
+   for (i=1; i<condZ; i++)
+     R->rayZ[i] = (cgcd(i,condZ) == 1)? isprincipalray(bnr, stoi(i)): gzero;
-     if (DEBUGLEVEL > 1 && prime[2] > cpt)
+   gptr[0] = &(R->L0);
-       { fprintferr("%ld ", prime[2]); cpt += 10; }
+   gptr[1] = &(R->L2);  gptr[2] = (GEN*)&(R->rayZ);
-   }
+   gptr[3] = &(R->L1);  gptr[5] = (GEN*)&(R->L1ray);
-   if (DEBUGLEVEL > 1) fprintferr("\n");
+   gptr[4] = &(R->L11); gptr[6] = (GEN*)&(R->L11ray);
+   gerepilemany(av,gptr,7);
+ }
-   for (i = 1; i <= cl; i++)
+ static void
-     CorrectCoeff((GEN)dataCR[i], matan[i], reduc[i], nmax, degs[i]);
+ InitPrimes(GEN bnr, long nmax, LISTray *R)
+ {
+   gpmem_t av = avma;
+   GEN bnf = (GEN)bnr[1], cond = gmael3(bnr,2,1,1);
+   long np,p,j, condZ = itos(gcoeff(cond,1,1));
+   GEN Npr, tabpr, prime, pr, nf = checknf(bnf);
+   byteptr d = diffptr + 1;
+   GEN *gptr[7];
-   FreeMat(matan2); FreeMat(reduc);
+   R->condZ = condZ;
-   avma = av; return matan;
+   R->L1 = cget1(nmax, t_VECSMALL);
+   R->L1ray = (GEN*)cget1(nmax, t_VEC);
+   prime = stoi(2);
+   for (p = 2; p <= nmax; prime[2] = p)
+   {
+     tabpr = primedec(nf, prime);
+     for (j = 1; j < lg(tabpr); j++)
+     {
+       pr  = (GEN)tabpr[j];
+       Npr = powgi((GEN)pr[1], (GEN)pr[4]);
+       if (is_bigint(Npr) || (np = Npr[2]) > nmax) continue;
+       if (condZ % p == 0 && idealval(nf, cond, pr)) continue;
+       appendL(R->L1, (GEN)np);
+       appendL((GEN)R->L1ray, isprincipalray(bnr, pr));
+     }
+     NEXT_PRIME_VIADIFF(p,d);
+   }
+   gptr[0] = &(R->L1);  gptr[1] = (GEN*)&(R->L1ray);
+   gerepilemany(av,gptr,2);
  }
- /********************************************************************/
- /*              5th part: compute L-functions at s=1                */
- /********************************************************************/
- /* if flag != 0, prec means decimal digits */
  static GEN
- get_limx(long r1, long r2, long prec, GEN *pteps, long flag)
+ get_limx(long r1, long r2, long prec, GEN *pteps)
  {
-   GEN eps, a, r, c0, A0, limx, Pi = mppi(prec), N, p1;
+   GEN eps, p1, a, c0, A0, limx, Pi2 = gmul2n(mppi(DEFAULTPREC), 1);
+   long r, N;
-   N = addss(r1, 2*r2);
+   N = r1 + 2*r2; r = r1 + r2;
-   a = gmul(gpow(gdeux, gsubgs(gdiv(stoi(r1), N), 1), DEFAULTPREC), N);
+   a = gmulgs(gpow(gdeux, gdiv(stoi(-2*r2), stoi(N)), DEFAULTPREC), N);
-   r = addss(r1, r2);
-   if (flag)
+   eps = real2n(-bit_accuracy(prec), DEFAULTPREC);
-     *pteps = eps = gmul2n(gpowgs(dbltor(10.), -prec), -1);
+   c0 = gpowgs(mpsqrt(Pi2), r-1);
-   else
+   c0 = gdivgs(gmul2n(c0,1), N);
-     *pteps = eps = gmul2n(gpowgs(dbltor(10.), (long)(-(prec-2) / pariK1)), -1);
+   c0 = gmul(c0, gpow(gdeux, gdiv(stoi(r1 * (r2-1)), stoi(2*N)),
-   c0 = gpow(gmul2n(Pi, 1), gdiv(subis(r, 1), gdeux), DEFAULTPREC);
-   c0 = gmul(c0, gdiv(gdeux, N));
-   c0 = gmul(c0, gpow(gdeux, gmul(gdiv(stoi(r1), gdeux),
-                                  gsubsg(1, gdiv(addis(r, 1), N))),
                       DEFAULTPREC));
-   A0 = glog(gdiv(gmul2n(c0, 1), eps), DEFAULTPREC);
+   A0 = glog(gdiv(gmul2n(c0,1), eps), DEFAULTPREC);
-   limx = gpow(gdiv(a, A0), gdiv(N, gdeux), DEFAULTPREC);
+   limx = gpow(gdiv(a, A0), gdiv(stoi(N), gdeux), DEFAULTPREC);
    p1   = gsub(glog(A0, DEFAULTPREC), glog(a, DEFAULTPREC));
-   p1   = gmul(gmul(p1, N), addis(r, 1));
+   p1   = gmulgs(p1, N*(r+1));
-   p1   = gdiv(p1, gmul2n(gadd(gmul2n(A0, 1), addis(r, 1)), 1));
+   p1   = gdiv(p1, gaddgs(gmul2n(A0, 2), 2*(r+1)));
-   limx = gmul(limx, gaddgs(p1, 1));
-   return limx;
+   if (pteps) *pteps = eps;
+   return gmul(limx, gaddgs(p1, 1));
  }
  static long
- GetBoundN0(GEN C,  long r1, long r2,  long prec, long flag)
+ GetBoundN0(GEN C,  long r1, long r2,  long prec)
  {
-   long av = avma;
+   gpmem_t av = avma;
-   GEN eps, limx = get_limx(r1, r2, prec, &eps, flag);
+   GEN limx = get_limx(r1, r2, prec, NULL);
    limx = gfloor(gdiv(C, limx));
    if (is_bigint(limx))
-Line 1724  GetBoundN0(GEN C,  long r1, long r2,  long prec, long
+Line 1642  GetBoundN0(GEN C,  long r1, long r2,  long prec, long
 Line 1724  GetBoundN0(GEN C,  long r1, long r2,  long prec, long
 Line 1642  GetBoundN0(GEN C,  long r1, long r2,  long prec, long
  static long
  GetBoundi0(long r1, long r2,  long prec)
  {
-   long av = avma, imin, i0, itest;
+   long imin, i0, itest;
-   GEN ftest, borneps, eps, limx = get_limx(r1, r2, prec, &eps, 0);
+   gpmem_t av = avma;
-   GEN Pi = mppi(DEFAULTPREC);
+   GEN ftest, borneps, eps, limx = get_limx(r1, r2, prec, &eps);
+   GEN sqrtPi = mpsqrt(mppi(DEFAULTPREC));
-   borneps = gmul(gmul2n(gun, r2), gpow(Pi, gdiv(subss(r2, 3), gdeux),
+   borneps = gmul(gmul2n(gun, r2), gpowgs(sqrtPi, r2-3));
-                                        DEFAULTPREC));
    borneps = gdiv(gmul(borneps, gpowgs(stoi(5), r1)), eps);
    borneps = gdiv(borneps, gsqrt(limx, DEFAULTPREC));
-Line 1740  GetBoundi0(long r1, long r2,  long prec)
+Line 1658  GetBoundi0(long r1, long r2,  long prec)
 Line 1740  GetBoundi0(long r1, long r2,  long prec)
 Line 1658  GetBoundi0(long r1, long r2,  long prec)
      itest = (i0 + imin) >> 1;
      ftest = gpowgs(limx, itest);
-     ftest = gmul(ftest, gpowgs(mpfactr(itest / 2, DEFAULTPREC), r1));
+     ftest = gmul(ftest, gpowgs(mpfactr(itest/2, DEFAULTPREC), r1));
-     ftest = gmul(ftest, gpowgs(mpfactr(itest, DEFAULTPREC), r2));
+     ftest = gmul(ftest, gpowgs(mpfactr(itest  , DEFAULTPREC), r2));
-     if(gcmp(ftest, borneps) >= 0)
+     if (gcmp(ftest, borneps) >= 0)
        i0 = itest;
      else
        imin = itest;
    }
    avma = av;
-   return (i0 / 2) * 2;
+   return i0 &= ~1; /* make it even */
  }
+ static GEN /* cf polcoeff */
+ _sercoeff(GEN x, long n)
+ {
+   long i = n - valp(x);
+   return (i < 0)? gzero: (GEN)x[i+2];
+ }
+ static void
+ affect_coeff(GEN q, long n, GEN y)
+ {
+   GEN x = _sercoeff(q,-n);
+   if (x == gzero) y[n] = zero; else gaffect(x, (GEN)y[n]);
+ }
+ typedef struct {
+   GEN c1, *aij, *bij, *powracpi, *cS, *cT;
+   long i0, a,b,c, r, rc1, rc2;
+ } ST_t;
  /* compute the principal part at the integers s = 0, -1, -2, ..., -i0
     of Gamma((s+1)/2)^a Gamma(s/2)^b Gamma(s)^c / (s - z) with z = 0 and 1 */
  /* NOTE: this is surely not the best way to do this, but it's fast enough! */
- static GEN
+ static void
- ppgamma(long a, long b, long c, long i0, long prec)
+ ppgamma(ST_t *T, long prec)
  {
-   GEN cst, gamun, gamdm, an, bn, cn_evn, cn_odd, x, x2, aij, p1, cf, p2;
+   GEN eul, gam,gamun,gamdm, an,bn,cn_evn,cn_odd, x,x2,X,Y, cf, sqpi;
-   long i, j, r, av = avma;
+   GEN p1, p2, *aij, *bij;
+   long a = T->a;
+   long b = T->b;
+   long c = T->c, r = T->r, i0 = T->i0;
+   long i,j, s,t;
+   gpmem_t av;
-   r = max(b + c + 1, a + c);
+   aij = (GEN*)cgetg(i0+1, t_VEC);
+   bij = (GEN*)cgetg(i0+1, t_VEC);
-   aij = cgetg(i0 + 1, t_VEC);
    for (i = 1; i <= i0; i++)
    {
-     aij[i] = lgetg(3, t_VEC);
+     aij[i] = p1 = cgetg(r+1, t_VEC);
-     mael(aij, i, 1) = lgetg(r + 1, t_VEC);
+     bij[i] = p2 = cgetg(r+1, t_VEC);
-     mael(aij, i, 2) = lgetg(r + 1, t_VEC);
+     for (j=1; j<=r; j++) { p1[j] = lgetr(prec); p2[j] = lgetr(prec); }
    }
+   av = avma;
    x   = polx[0];
-   x2  = gmul2n(x, -1);
+   x2  = gmul2n(x, -1); /* x/2 */
+   eul = mpeuler(prec);
+   sqpi= mpsqrt(mppi(prec)); /* Gamma(1/2) */
-   /* Euler gamma constant, values of Riemann zeta functions at
+   /* expansion of log(Gamma(u)) at u = 1 */
-      positive integers */
+   gamun = cgetg(r+3, t_SER);
-   cst = cgetg(r + 2, t_VEC);
-   cst[1] = (long)mpeuler(prec);
-   for (i = 2; i <= r + 1; i++)
-     cst[i] = (long)gzeta(stoi(i), prec);
-   /* the expansion of log(Gamma(s)) at s = 1 */
-   gamun = cgetg(r + 2, t_SER);
    gamun[1] = evalsigne(1) | evalvalp(0) | evalvarn(0);
    gamun[2] = zero;
-   for (i = 1; i <= r; i++)
+   gamun[3] = lneg(eul);
+   for (i = 2; i <= r; i++)
    {
-     gamun[i + 2] = ldivgs((GEN)cst[i], i);
+     p1 = gdivgs(gzeta(stoi(i),prec), i);
-     if (i%2) gamun[i + 2] = lneg((GEN)gamun[i + 2]);
+     if (odd(i)) p1 = gneg(p1);
+     gamun[i+2] = (long)p1;
    }
+   gamun = gexp(gamun, prec); /* Gamma(1 + x) */
+   gam = gdiv(gamun,x); /* Gamma(x) */
-   /* the expansion of log(Gamma(s)) at s = 1/2 */
+   /* expansion of log(Gamma(u) / Gamma(1/2)) at u = 1/2 */
-   gamdm = cgetg(r + 2, t_SER);
+   gamdm = cgetg(r+3, t_SER);
    gamdm[1] = evalsigne(1) | evalvalp(0) | evalvarn(0);
-   gamdm[2] = (long)mplog(gsqrt(mppi(prec), prec));
+   gamdm[2] = zero;
-   gamdm[3] = lneg(gadd(gmul2n(glog(gdeux, prec), 1), (GEN)cst[1]));
+   gamdm[3] = lneg(gadd(gmul2n(mplog2(prec), 1), eul));
    for (i = 2; i <= r; i++)
-     gamdm[i + 2] = lmul((GEN)gamun[i + 2], subis(shifti(gun, i), 1));
+     gamdm[i+2] = lmul((GEN)gamun[i+2], subis(shifti(gun,i), 1));
+   gamdm = gmul(sqpi, gexp(gamdm, prec)); /* Gamma(1/2 + x) */
-   gamun = gexp(gamun, prec);
+  /* We simplify to get one of the following two expressions
-   gamdm = gexp(gamdm, prec);
+   * if (b > a) : sqrt{Pi}^a 2^{a-au} Gamma(u)^{a+c} Gamma(  u/2  )^{|b-a|}
+   * if (b <= a): sqrt{Pi}^b 2^{b-bu} Gamma(u)^{b+c) Gamma((u+1)/2)^{|b-a|} */
-   /* We simplify to get one of the following two expressions */
-   /* Case 1 (b > a): sqrt{Pi}^a 2^{a - as} Gamma(s/2)^{b-a} Gamma(s)^{a + c} */
    if (b > a)
    {
-     cf = gpui(mppi(prec), gmul2n(stoi(a), -1), prec);
+     t = a; s = b; X = x2; Y = gsub(x2,ghalf);
+     p1 = gsubst(gam,0,x2);
-     /* an is the expansion of Gamma(x)^{a+c} */
+     p2 = gdiv(gsubst(gamdm,0,x2), Y); /* Gamma((x-1)/2) */
-     an = gpowgs(gdiv(gamun, x), a + c);
-     /* bn is the expansion of 2^{a - ax} */
-     bn = gpowgs(gpow(gdeux, gsubsg(1, x), prec), a);
-     /* cn_evn is the expansion of Gamma(x/2)^{b-a} */
-     cn_evn = gpowgs(gdiv(gsubst(gamun, 0, x2), x2), b - a);
-     /* cn_odd is the expansion of Gamma((x-1)/2)^{b-a} */
-     cn_odd = gpowgs(gdiv(gsubst(gamdm, 0, x2), gsub(x2, ghalf)), b - a);
-     for (i = 0; i < i0/2; i++)
-     {
-       p1 = gmul(cf, gmul(an, gmul(bn, cn_evn)));
-       p2 = gdiv(p1, gsubgs(x, 2*i));
-       for (j = 1; j <= r; j++)
-         mael3(aij, 2*i + 1, 1, j) = (long)polcoeff0(p2, -j, 0);
-       p2 = gdiv(p1, gsubgs(x, 2*i + 1));
-       for (j = 1; j <= r; j++)
-         mael3(aij, 2*i + 1, 2, j) = (long)polcoeff0(p2, -j, 0);
-       /* an(x-s-1) = an(x-s) / (x-s-1)^{a+c} */
-       an = gdiv(an, gpowgs(gsubgs(x, 2*i + 1), a + c));
-       /* bn(x-s-1) = 2^a bn(x-s) */
-       bn = gmul2n(bn, a);
-       /* cn_evn(x-s-2) = cn_evn(x-s) / (x/2 - (s+2)/2)^{b-a} */
-       cn_evn = gdiv(cn_evn, gpowgs(gsubgs(x2, i + 1), b - a));
-       p1 = gmul(cf, gmul(an, gmul(bn, cn_odd)));
-       p2 = gdiv(p1, gsubgs(x, 2*i + 1));
-       for (j = 1; j <= r; j++)
-         mael3(aij, 2*i + 2, 1, j) = (long)polcoeff0(p2, -j, 0);
-       p2 = gdiv(p1, gsubgs(x, 2*i + 2));
-       for (j = 1; j <= r; j++)
-         mael3(aij, 2*i + 2, 2, j) = (long)polcoeff0(p2, -j, 0);
-       an = gdiv(an, gpowgs(gsubgs(x, 2*i + 2), a + c));
-       bn = gmul2n(bn, a);
-       /* cn_odd(x-s-2) = cn_odd(x-s) / (x/2 - (s+2)/2)^{b-a} */
-       cn_odd = gdiv(cn_odd, gpowgs(gsub(x2, gaddgs(ghalf, i + 1)), b - a));
-     }
    }
    else
-   /* Case 2 (b <= a): sqrt{Pi}^b 2^{b - bs} Gamma((s+1)/2)^{a-b}
-                                                          Gamma(s)^{b + c) */
    {
-     cf = gpui(mppi(prec), gmul2n(stoi(b), -1), prec);
+     t = b; s = a; X = gadd(x2,ghalf); Y = x2;
+     p1 = gsubst(gamdm,0,x2);
+     p2 = gsubst(gam,0,x2);
+   }
+   cf = gpowgs(sqpi, t);
+   an = gpowgs(gpow(gdeux, gsubsg(1,x), prec), t); /* 2^{t-tx} */
+   bn = gpowgs(gam, t+c); /* Gamma(x)^{t+c} */
+   cn_evn = gpowgs(p1, s-t); /* Gamma(X)^{s-t} */
+   cn_odd = gpowgs(p2, s-t); /* Gamma(Y)^{s-t} */
+   for (i = 0; i < i0/2; i++)
+   {
+     GEN C1,q1, A1 = aij[2*i+1], B1 = bij[2*i+1];
+     GEN C2,q2, A2 = aij[2*i+2], B2 = bij[2*i+2];
-     /* an is the expansion of Gamma(x)^{b+c} */
+     C1 = gmul(cf, gmul(bn, gmul(an, cn_evn)));
-     an = gpowgs(gdiv(gamun, x), b + c);
+     p1 = gdiv(C1, gsubgs(x, 2*i));
+     q1 = gdiv(C1, gsubgs(x, 2*i+1));
-     /* bn is the expansion of 2^{b - bx} */
+     /* an(x-u-1) = 2^t an(x-u) */
-     bn = gpowgs(gpow(gdeux, gsubsg(1, x), prec), b);
+     an = gmul2n(an, t);
+     /* bn(x-u-1) = bn(x-u) / (x-u-1)^{t+c} */
+     bn = gdiv(bn, gpowgs(gsubgs(x, 2*i+1), t+c));
-     /* cn_evn is the expansion of Gamma((x+1)/2)^{a-b} */
+     C2 = gmul(cf, gmul(bn, gmul(an, cn_odd)));
-     cn_evn = gpowgs(gsubst(gamdm, 0, x2), a - b);
+     p2 = gdiv(C2, gsubgs(x, 2*i+1));
+     q2 = gdiv(C2, gsubgs(x, 2*i+2));
-     /* cn_odd is the expansion of Gamma(x/2)^{a-b} */
+     for (j = 1; j <= r; j++)
-     cn_odd = gpowgs(gdiv(gsubst(gamun, 0, x2), x2), a - b);
-     for (i = 0; i < i0/2; i++)
      {
-       p1 = gmul(cf, gmul(an, gmul(bn, cn_evn)));
+       affect_coeff(p1, j, A1); affect_coeff(q1, j, B1);
+       affect_coeff(p2, j, A2); affect_coeff(q2, j, B2);
+     }
-       p2 = gdiv(p1, gsubgs(x, 2*i));
+     an = gmul2n(an, t);
-       for (j = 1; j <= r; j++)
+     bn = gdiv(bn, gpowgs(gsubgs(x, 2*i+2), t+c));
-         mael3(aij, 2*i + 1, 1, j) = (long)polcoeff0(p2, -j, 0);
+     /* cn_evn(x-2i-2) = cn_evn(x-2i)  / (X - (i+1))^{s-t} */
+     /* cn_odd(x-2i-3) = cn_odd(x-2i-1)/ (Y - (i+1))^{s-t} */
+     cn_evn = gdiv(cn_evn, gpowgs(gsubgs(X,i+1), s-t));
+     cn_odd = gdiv(cn_odd, gpowgs(gsubgs(Y,i+1), s-t));
+   }
+   T->aij = aij;
+   T->bij = bij; avma = av;
+ }
-       p2 = gdiv(p1, gsubgs(x, 2*i + 1));
+ static GEN
-       for (j = 1; j <= r; j++)
+ _cond(GEN dtcr, int quad)
-         mael3(aij, 2*i + 1, 2, j) = (long)polcoeff0(p2, -j, 0);
+ {
+   GEN cond;
+   if (quad) cond = (GEN)dtcr[7];
+   else
+   {
+     cond = cgetg(3, t_VEC);
+     cond[1] = dtcr[7];
+     cond[2] = dtcr[9];
+   }
+   return cond;
+ }
-       /* an(x-s-1) = an(x-s) / (x-s-1)^{b+c} */
+ /* sort chars according to conductor */
-       an = gdiv(an, gpowgs(gsubgs(x, 2*i + 1), b + c));
+ static GEN
+ sortChars(GEN dataCR, int quad)
+ {
+   const long cl = lg(dataCR) - 1;
+   GEN vCond  = cgetg(cl+1, t_VEC);
+   GEN CC     = cgetg(cl+1, t_VECSMALL);
+   GEN nvCond = cgetg(cl+1, t_VECSMALL);
+   long j,k, ncond;
+   GEN vChar;
-       /* bn(x-s-1) = 2^b bn(x-s) */
+   for (j = 1; j <= cl; j++) nvCond[j] = 0;
-       bn = gmul2n(bn, b);
-       /* cn_evn(x-s-2) = cn_evn(x-s) / (x/2 - (s+1)/2)^{a-b} */
+   ncond = 0;
-       cn_evn = gdiv(cn_evn, gpowgs(gsub(x2, gaddgs(ghalf, i)), a - b));
+   for (j = 1; j <= cl; j++)
+   {
+     GEN cond = _cond((GEN)dataCR[j], quad);
+     for (k = 1; k <= ncond; k++)
+       if (gegal(cond, (GEN)vCond[k])) break;
+     if (k > ncond) vCond[++ncond] = (long)cond;
+     nvCond[k]++; CC[j] = k; /* char j has conductor number k */
+   }
+   vChar = cgetg(ncond+1, t_VEC);
+   for (k = 1; k <= ncond; k++)
+   {
+     vChar[k] = lgetg(nvCond[k]+1, t_VECSMALL);
+     nvCond[k] = 0;
+   }
+   for (j = 1; j <= cl; j++)
+   {
+     k = CC[j]; nvCond[k]++;
+     mael(vChar,k,nvCond[k]) = j;
+   }
+   return vChar;
+ }
-       p1 = gmul(cf, gmul(an, gmul(bn, cn_odd)));
+ static void
+ get_cS_cT(ST_t *T, long n)
+ {
+   gpmem_t av;
+   GEN csurn, nsurc, lncsurn;
+   GEN A,B,s,t,Z,*aij,*bij;
+   long i,j,r,i0;
-       p2 = gdiv(p1, gsubgs(x, 2*i + 1));
+   if (T->cS[n]) return;
-       for (j = 1; j <= r; j++)
-         mael3(aij, 2*i + 2, 1, j) = (long)polcoeff0(p2, -j, 0);
-       p2 = gdiv(p1, gsubgs(x, 2*i + 2));
+   av = avma;
-       for (j = 1; j <= r; j++)
+   aij = T->aij; i0= T->i0;
-         mael3(aij, 2*i + 2, 2, j) = (long)polcoeff0(p2, -j, 0);
+   bij = T->bij; r = T->r;
+   Z = cgetg(r+1, t_VEC);
+   Z[1] = un;
-       an = gdiv(an, gpowgs(gsubgs(x, 2*i + 2), b + c));
+   csurn = gdivgs(T->c1, n);
-       bn = gmul2n(bn, b);
+   nsurc = ginv(csurn);
+   lncsurn = mplog(csurn);
-       /* cn_odd(x-s-2) = cn_odd(x-s) / (x/2 - (s+1)/2)^{a-b} */
+   Z[2] = (long)lncsurn; /* r >= 2 */
-       cn_odd = gdiv(cn_odd, gpowgs(gsubgs(x2, i + 1), a - b));
+   for (i = 3; i <= r; i++)
+   {
+     s = gmul((GEN)Z[i-1], lncsurn);
+     Z[i] = ldivgs(s, i-1); /* Z[i] = ln^(i-1)(c1/n) / (i-1)! */
+   }
+   /* i = i0 */
+     A = aij[i0]; t = (GEN)A[1];
+     B = bij[i0]; s = (GEN)B[1];
+     for (j = 2; j <= r; j++)
+     {
+       s = gadd(s, gmul((GEN)Z[j], (GEN)B[j]));
+       t = gadd(t, gmul((GEN)Z[j], (GEN)A[j]));
      }
+   for (i = i0 - 1; i > 1; i--)
+   {
+     A = aij[i]; t = gmul(t, nsurc);
+     B = bij[i]; s = gmul(s, nsurc);
+     for (j = odd(i)? T->rc2: T->rc1; j; j--)
+     {
+       s = gadd(s, gmul((GEN)Z[j], (GEN)B[j]));
+       t = gadd(t, gmul((GEN)Z[j], (GEN)A[j]));
+     }
    }
+   /* i = 1 */
+     A = aij[1]; t = gmul(t, nsurc);
+     B = bij[1]; s = gmul(s, nsurc);
+     for (j = 1; j <= r; j++)
+     {
+       s = gadd(s, gmul((GEN)Z[j], (GEN)B[j]));
+       t = gadd(t, gmul((GEN)Z[j], (GEN)A[j]));
+     }
+   s = gadd(s, gmul(csurn, T->powracpi[T->b]));
+   T->cS[n] = gclone(s);
+   T->cT[n] = gclone(t); avma = av;
+ }
-   return gerepilecopy(av, aij);
+ static void
+ clear_cScT(ST_t *T, long N)
+ {
+   GEN *cS = T->cS, *cT = T->cT;
+   long i;
+   for (i=1; i<=N; i++)
+     if (cS[i]) { gunclone(cS[i]); gunclone(cT[i]); cS[i] = cT[i] = NULL; }
  }
+ static void
+ init_cScT(ST_t *T, GEN CHI, long N, long prec)
+ {
+   GEN p1 = (GEN)CHI[9];
+   T->a = p1[1];
+   T->b = p1[2];
+   T->c = p1[3];
+   T->rc1 = T->a + T->c;
+   T->rc2 = T->b + T->c;
+   T->r   = max(T->rc2+1, T->rc1); /* >= 2 */
+   ppgamma(T, prec);
+   clear_cScT(T, N);
+ }
  static GEN
- GetST(GEN dataCR, long prec)
+ GetST(GEN dataCR, GEN vChar, long prec)
  {
-   GEN N0, CC, bnr, bnf, nf, Pi, racpi, C, cond, aij, B, S, T, csurn, lncsurn;
+   const long cl = lg(dataCR) - 1;
-   GEN degs, p1, p2, nsurc, an, rep, powlncn, powracpi;
+   gpmem_t av, av1, av2;
-   long i, j, k, n, av = avma, av1, av2, hk, fj, id, prec2, i0, nmax;
+   long ncond, n, j, k, jc, nmax, prec2, i0, r1, r2;
-   long a, b, c, rc1, rc2, r, r1, r2;
+   GEN bnr, nf, racpi, *powracpi;
-   int ***matan;
+   GEN rep, N0, C, T, S, an, degs;
+   LISTray LIST;
+   ST_t cScT;
-   if (DEBUGLEVEL) timer2();
+   bnr = gmael(dataCR,1,4);
-   bnr   = gmael(dataCR, 1, 4);
+   nf  = checknf(bnr);
-   bnf   = checkbnf(bnr);
+   /* compute S,T differently if nf is quadratic */
-   nf    = checknf(bnf);
+   if (degpol(nf[1]) == 2) return QuadGetST(dataCR, vChar, prec);
-   r1    = nf_get_r1(nf);
-   r2    = nf_get_r2(nf);
-   hk    = lg(dataCR) - 1;
-   prec2 = ((prec - 2)<<1) + EXTRA_PREC;
-   Pi    = mppi(prec2);
+   if (DEBUGLEVEL) (void)timer2();
-   racpi = gsqrt(Pi, prec2);
+   /* allocate memory for answer */
+   rep = cgetg(3, t_VEC);
+   S = cgetg(cl+1, t_VEC); rep[1] = (long)S;
+   T = cgetg(cl+1, t_VEC); rep[2] = (long)T;
+   for (j = 1; j <= cl; j++)
+   {
+     S[j] = (long)cgetc(prec);
+     T[j] = (long)cgetc(prec);
+   }
+   av = avma;
-   C    = cgetg(hk + 1, t_VEC);
+   /* initializations */
-   cond = cgetg(hk + 1, t_VEC);
+   degs = GetDeg(dataCR);
-   N0 = new_chunk(hk+1);
+   ncond = lg(vChar)-1;
-   CC = new_chunk(hk+1);
+   nf_get_sign(nf,&r1,&r2);
+   C  = cgetg(ncond+1, t_VEC);
+   N0 = cgetg(ncond+1, t_VECSMALL);
    nmax = 0;
-   for (i = 1; i <= hk; i++)
+   for (j = 1; j <= ncond; j++)
    {
-     C[i]    = mael(dataCR, i, 2);
+     C[j]  = mael(dataCR, mael(vChar,j,1), 2);
+     N0[j] = GetBoundN0((GEN)C[j], r1, r2, prec);
-     p1 = cgetg(3, t_VEC);
+     if (nmax < N0[j]) nmax  = N0[j];
-     p1[1] = mael(dataCR, i, 7);
-     p1[2] = mael(dataCR, i, 9);
-     cond[i] = (long)p1;
-     N0[i] = GetBoundN0((GEN)C[i], r1, r2, prec, 0);
-     if (nmax < N0[i]) nmax  = N0[i];
    }
    if ((ulong)nmax > maxprime())
-     err(talker, "Not enough precomputed primes (need all primes up to %ld)", nmax);
+     err(talker, "Not enough precomputed primes (need all p <= %ld)", nmax);
    i0 = GetBoundi0(r1, r2, prec);
-   if(DEBUGLEVEL > 1) fprintferr("nmax = %ld and i0 = %ld\n", nmax, i0);
+   if (DEBUGLEVEL>1) fprintferr("nmax = %ld, i0 = %ld\n", nmax, i0);
+   InitPrimes(gmael(dataCR,1,4), nmax, &LIST);
-   matan = ComputeCoeff(dataCR, nmax, prec);
+   prec2 = ((prec-2) << 1) + EXTRA_PREC;
-   degs = GetDeg(dataCR);
+   racpi = mpsqrt(mppi(prec2));
-   if (DEBUGLEVEL) msgtimer("Compute an");
+   powracpi = (GEN*)cgetg(r1+2,t_VEC);
+   powracpi++; powracpi[0] = gun; powracpi[1] = racpi;
+   for (j=2; j<=r1; j++) powracpi[j] = mulrr(powracpi[j-1], racpi);
+   cScT.powracpi = powracpi;
-   p1 = cgetg(3, t_COMPLEX);
+   cScT.cS = (GEN*)cgetg(nmax+1, t_VEC);
-   p1[1] = lgetr(prec2);
+   cScT.cT = (GEN*)cgetg(nmax+1, t_VEC);
-   p1[2] = lgetr(prec2);
+   for (j=1; j<=nmax; j++) cScT.cS[j] = cScT.cT[j] = NULL;
-   gaffect(gzero, p1);
-   S = cgetg(hk + 1, t_VEC);
+   cScT.i0 = i0;
-   T = cgetg(hk + 1, t_VEC);
-   for (id = 1; id <= hk; id++)
+   av1 = avma;
+   for (jc = 1; jc <= ncond; jc++)
    {
-     S[id] = lcopy(p1);
+     const GEN LChar = (GEN)vChar[jc];
-     T[id] = lcopy(p1);
+     const long nChar = lg(LChar)-1, NN = N0[jc];
-     for (k = 1; k < id; k++)
-       if (gegal((GEN)cond[id], (GEN)cond[k])) break;
-     CC[id] = k;
-   }
-   powracpi = cgetg(hk + 1, t_VEC);
+     if (DEBUGLEVEL>1)
-   for (j = 1; j <= hk; j++)
+       fprintferr("* conductor no %ld/%ld (N = %ld)\n\tInit: ", jc,ncond,NN);
-     powracpi[j] = (long)gpow(racpi, gmael3(dataCR, j, 9, 2), prec2);
-   av1 = avma;
+     cScT.c1 = (GEN)C[jc];
-   if (DEBUGLEVEL > 1) fprintferr("n = ");
+     init_cScT(&cScT, (GEN)dataCR[LChar[1]], NN, prec2);
-   for (id = 1; id <= hk; id++)
-   {
-     if (CC[id] != id) continue;
-     p2 = gmael(dataCR, id, 9);
-     a  = itos((GEN)p2[1]);
-     b  = itos((GEN)p2[2]);
-     c  = itos((GEN)p2[3]);
-     aij = ppgamma(a, b, c, i0, prec2);
-     rc1 = a + c;
-     rc2 = b + c; r = max(rc2 + 1, rc1);
      av2 = avma;
+     for (k = 1; k <= nChar; k++)
-     for (n = 1; n <= N0[id]; n++)
      {
-       if (DEBUGLEVEL > 1 && n%100 == 0) fprintferr("%ld ", n);
+       const long t = LChar[k], d = degs[t];
+       const GEN z = gmael3(dataCR, t, 5, 2);
-       for (k = 1; k <= hk; k++)
+       GEN p1 = gzero, p2 = gzero;
-         if (CC[k] == id && !IsZero(matan[k][n], degs[k])) break;
+       long c = 0;
-       if (k > hk) continue;
+       int **matan;
-       csurn = gdivgs((GEN)C[id], n);
+       if (DEBUGLEVEL>1)
-       nsurc = ginv(csurn);
+         fprintferr("\tcharacter no: %ld (%ld/%ld)\n", t,k,nChar);
+       matan = ComputeCoeff((GEN)dataCR[t], &LIST, NN, d);
-       B = cgetg(r + 1, t_VEC);
+       for (n = 1; n <= NN; n++)
-       lncsurn = glog(csurn, prec2);
+         if ((an = EvalCoeff(z, matan[n], d)))
-       powlncn = gun;
-       fj = 1;
-       p1 = gzero;
-       p2 = gzero;
-       for (j = 1; j <= r; j++)
-       {
-         if (j > 2) fj = fj * (j - 1);
-         B[j] = ldivgs(powlncn, fj);
-         p1 = gadd(p1, gmul((GEN)B[j], gmael3(aij, i0, 2, j)));
-         p2 = gadd(p2, gmul((GEN)B[j], gmael3(aij, i0, 1, j)));
-         powlncn = gmul(powlncn, lncsurn);
-       }
-       for (i = i0 - 1; i > 1; i--)
-       {
-         p1 = gmul(p1, nsurc);
-         p2 = gmul(p2, nsurc);
-         for (j = i%2? rc2: rc1; j; j--)
          {
-           p1 = gadd(p1, gmul((GEN)B[j], gmael3(aij, i, 2, j)));
+           get_cS_cT(&cScT, n);
-           p2 = gadd(p2, gmul((GEN)B[j], gmael3(aij, i, 1, j)));
+           p1 = gadd(p1, gmul(an,        cScT.cS[n]));
+           p2 = gadd(p2, gmul(gconj(an), cScT.cT[n]));
+           if (++c == 256)
+           { GEN *gptr[2]; gptr[0]=&p1; gptr[1]=&p2;
+             gerepilemany(av2,gptr,2); c = 0;
+           }
          }
-       }
+       gaffect(p1, (GEN)S[t]);
-       p1 = gmul(p1, nsurc);
+       gaffect(p2, (GEN)T[t]);
-       p2 = gmul(p2, nsurc);
+       FreeMat(matan, NN); avma = av2;
-       for (j = 1; j <= r; j++)
-       {
-         p1 = gadd(p1, gmul((GEN)B[j], gmael3(aij, 1, 2, j)));
-         p2 = gadd(p2, gmul((GEN)B[j], gmael3(aij, 1, 1, j)));
-       }
-       p1 = gadd(p1, gmul(csurn, (GEN)powracpi[id]));
-       for (j = 1; j <= hk; j++)
-         if (CC[j] == id &&
-             (an = EvalCoeff(gmael3(dataCR, j, 5, 2), matan[j][n], degs[j])))
-         {
-           gaffect(gadd((GEN)S[j], gmul(p1, an)),        (GEN)S[j]);
-           gaffect(gadd((GEN)T[j], gmul(p2, gconj(an))), (GEN)T[j]);
-         }
-       avma = av2;
      }
+     if (DEBUGLEVEL>1) fprintferr("\n");
      avma = av1;
    }
-   FreeMat(matan);
+   if (DEBUGLEVEL) msgtimer("S&T");
+   clear_cScT(&cScT, nmax);
-   if (DEBUGLEVEL > 1) fprintferr("\n");
+   avma = av; return rep;
-   if (DEBUGLEVEL) msgtimer("Compute S&T");
-   rep = cgetg(3, t_VEC);
-   rep[1] = (long)S;
-   rep[2] = (long)T;
-   return gerepilecopy(av, rep);
  }
- /* Given datachi, S(chi) and T(chi), return L(1, chi) if fl = 1,
+ /* Given W(chi) [possibly NULL], S(chi) and T(chi), return L(1, chi) if fl = 1,
     or [r(chi), c(chi)] where r(chi) is the rank of chi and c(chi)
     is given by L(s, chi) = c(chi).s^r(chi) at s = 0 if fl = 0.
     if fl2 = 1, adjust the value to get L_S(s, chi). */
  static GEN
- GetValue(GEN datachi, GEN S, GEN T, long fl, long fl2, long prec)
+ GetValue(GEN dtcr, GEN W, GEN S, GEN T, long fl, long fl2, long prec)
  {
-   GEN W, A, q, b, c, d, rchi, cf, VL, rep, racpi, nS, nT;
+   gpmem_t av = avma;
-   long av = avma;
+   GEN A, cf, VL, rep, racpi, p1;
+   long q, b, c;
+   int isreal;
-   racpi = gsqrt(mppi(prec), prec);
+   racpi = mpsqrt(mppi(prec));
-   W = ComputeArtinNumber(datachi, 0, prec);
+   if (!W)
-   A = ComputeAChi(datachi, fl, prec);
+     W = (GEN)ComputeArtinNumber((GEN)dtcr[3], _vec((GEN)dtcr[8]), 1, prec)[1];
-   d = gmael(datachi, 8, 3);
+   isreal = (itos(gmael(dtcr, 8, 3)) <= 2);
-   q = gmael(datachi, 9, 1);
+   p1 = (GEN)dtcr[9];
-   b = gmael(datachi, 9, 2);
+   q = p1[1];
-   c = gmael(datachi, 9, 3);
+   b = p1[2];
+   c = p1[3];
-   rchi = addii(b, c);
    if (!fl)
-   {
+   { /* VL = (W(chi).S(conj(chi)) + T(chi)) / (sqrt(Pi)^q 2^{r1 - q}) */
-     cf = gmul2n(gpow(racpi, q, 0), itos(b));
+     GEN rchi = stoi(b + c);
+     cf = gmul2n(gpowgs(racpi, q), b);
-     nS = gdiv(gconj(S), cf);
+     VL = gadd(gmul(W, gconj(S)), gconj(T));
-     nT = gdiv(gconj(T), cf);
+     if (isreal) VL = greal(VL);
+     VL = gdiv(VL, cf);
-     /* VL = W(chi).S(conj(chi)) + T(chi)) / (sqrt(Pi)^q 2^{r1 - q}) */
-     VL = gadd(gmul(W, nS), nT);
-     if (cmpis(d, 3) < 0) VL = greal(VL);
      if (fl2)
      {
+       A = ComputeAChi(dtcr, fl, prec);
        VL = gmul((GEN)A[2], VL);
        rchi = gadd(rchi, (GEN)A[1]);
      }
-Line 2119  GetValue(GEN datachi, GEN S, GEN T, long fl, long fl2,
+Line 2079  GetValue(GEN datachi, GEN S, GEN T, long fl, long fl2,
 Line 2119  GetValue(GEN datachi, GEN S, GEN T, long fl, long fl2,
 Line 2079  GetValue(GEN datachi, GEN S, GEN T, long fl, long fl2,
      rep[2] = (long)VL;
    }
    else
-   {
+   { /* VL = S(chi) + W(chi).T(chi)) / (C(chi) sqrt(Pi)^{r1 - q}) */
-     cf = gmul((GEN)datachi[2], gpow(racpi, b, 0));
+     cf = gmul((GEN)dtcr[2], gpowgs(racpi, b));
-     /* VL = S(chi) + W(chi).T(chi)) / (C(chi) sqrt(Pi)^{r1 - q}) */
+     rep = gadd(S, gmul(W, T));
-     rep = gdiv(gadd(S, gmul(W, T)), cf);
+     if (isreal) rep = greal(rep);
-     if (cmpis(d, 3) < 0) rep = greal(rep);
+     rep = gdiv(rep, cf);
-     if (fl2) rep = gmul(A, rep);
+     if (fl2)
+     {
+       A = ComputeAChi(dtcr, fl, prec);
+       rep = gmul(A, rep);
+     }
    }
-  return gerepilecopy(av, rep);
+   return gerepilecopy(av, rep);
  }
  /* return the order and the first non-zero term of L(s, chi0)
-Line 2139  GetValue1(GEN bnr, long flag, long prec)
+Line 2103  GetValue1(GEN bnr, long flag, long prec)
 Line 2139  GetValue1(GEN bnr, long flag, long prec)
 Line 2103  GetValue1(GEN bnr, long flag, long prec)
  {
    GEN bnf = checkbnf(bnr), nf = checknf(bnf);
    GEN hk, Rk, wk, c, rep, mod0, diff;
-   long i, l, r, r1, r2, av = avma;
+   long i, l, r, r1, r2;
+   gpmem_t av = avma;
    r1 = nf_get_r1(nf);
    r2 = nf_get_r2(nf);
    hk = gmael3(bnf, 8, 1, 1);
    Rk = gmael(bnf, 8, 2);
    wk = gmael3(bnf, 8, 4, 1);
    c = gneg_i(gdiv(gmul(hk, Rk), wk));
    r = r1 + r2 - 1;
-Line 2163  GetValue1(GEN bnr, long flag, long prec)
+Line 2128  GetValue1(GEN bnr, long flag, long prec)
 Line 2163  GetValue1(GEN bnr, long flag, long prec)
 Line 2128  GetValue1(GEN bnr, long flag, long prec)
    rep = cgetg(3, t_VEC);
    rep[1] = lstoi(r);
    rep[2] = (long)c;
    return gerepilecopy(av, rep);
  }
  /********************************************************************/
  /*                6th part: recover the coefficients                */
  /********************************************************************/
  static long
- TestOne(GEN plg,  GEN beta,  GEN B,  long v,  long G,  long N)
+ TestOne(GEN plg, RC_data *d)
  {
-   long j;
+   long j, v = d->v;
    GEN p1;
-   p1 = gsub(beta, (GEN)plg[v]);
+   p1 = gsub(d->beta, (GEN)plg[v]);
-   if (expo(p1) >= G) return 0;
+   if (expo(p1) >= d->G) return 0;
-   for (j = 1; j <= N; j++)
+   for (j = 1; j <= d->N; j++)
      if (j != v)
      {
        p1 = gabs((GEN)plg[j], DEFAULTPREC);
-       if (gcmp(p1, B) > 0) return 0;
+       if (gcmp(p1, d->B) > 0) return 0;
      }
    return 1;
  }
- /* Using linear dependance relations */
  static GEN
- RecCoeff2(GEN nf,  GEN beta,  GEN B,  long v,  long prec)
+ chk_reccoeff_init(FP_chk_fun *chk, GEN gram, GEN mat)
  {
-   long N, G, i, bacmin, bacmax, av = avma, av2;
+   RC_data *d = (RC_data*)chk->data;
-   GEN vec, velt, p1, cand, M, plg, pol, cand2;
+   (void)gram; d->U = mat; return d->nB;
-   M    = gmael(nf, 5, 1);
-   pol  = (GEN)nf[1];
-   N    = degpol(pol);
-   vec  = gtrans((GEN)gtrans(M)[v]);
-   velt = (GEN)nf[7];
-   G = min( - 20, - bit_accuracy(prec) >> 4);
-   p1 = cgetg(2, t_VEC);
-   p1[1] = lneg(beta);
-   vec = concat(p1, vec);
-   p1[1] = zero;
-   velt = concat(p1, velt);
-   bacmin = (long)(.225 * bit_accuracy(prec));
-   bacmax = (long)(.315 * bit_accuracy(prec));
-   av2 = avma;
-   for (i = bacmax; i >= bacmin; i--)
-   {
-     p1 = lindep2(vec, i);
-     if (signe((GEN)p1[1]))
-     {
-       p1    = ground(gdiv(p1, (GEN)p1[1]));
-       cand  = gmodulcp(gmul(velt, gtrans(p1)), pol);
-       cand2 = algtobasis(nf, cand);
-       plg   = gmul(M, cand2);
-       if (TestOne(plg, beta, B, v, G, N))
-         return gerepilecopy(av, cand);
-     }
-     avma = av2;
-   }
-   return NULL;
  }
- GEN
+ static GEN
- chk_reccoeff_init(FP_chk_fun *chk, GEN nf, GEN gram, GEN mat, long *ptprec)
+ chk_reccoeff(void *data, GEN x)
  {
-   GEN data = chk->data;
+   RC_data *d = (RC_data*)data;
-   data[6] = (long)mat;
+   long N = d->N, j;
-   chk->data = data;
+   GEN p1 = gmul(d->U, x), sol, plg;
-   return (GEN)data[7];
- }
- GEN
- chk_reccoeff(GEN data, GEN x)
- {
-   GEN M = (GEN)data[0], beta = (GEN)data[1], B = (GEN)data[2];
-   long v = data[3], G = data[4], N = data[5], j;
-   GEN U = (GEN)data[6], p1 = gmul(U, x), sol, plg;
    if (!gcmp1((GEN)p1[1])) return NULL;
    sol = cgetg(N + 1, t_COL);
    for (j = 1; j <= N; j++)
      sol[j] = lmulii((GEN)p1[1], (GEN)p1[j + 1]);
-   plg = gmul(M, sol);
+   plg = gmul(d->M, sol);
-   if (TestOne(plg, beta, B, v, G, N)) return sol;
+   if (TestOne(plg, d)) return sol;
    return NULL;
  }
- GEN
- chk_reccoeff_post(GEN data, GEN res)
- {
-   return res;
- }
  /* Using Cohen's method */
  static GEN
- RecCoeff3(GEN nf, GEN beta, GEN B, long v, long prec)
+ RecCoeff3(GEN nf, RC_data *d, long prec)
  {
-   GEN A, M, nB, cand, p1, B2, C2, data, tB, beta2, eps, nf2, Bd;
+   GEN A, M, nB, cand, p1, B2, C2, tB, beta2, eps, nf2, Bd;
-   long N, G, i, j, k, l, ct = 0, av = avma, prec2;
+   GEN beta = d->beta, B = d->B;
-   FP_chk_fun *chk;
+   long N = d->N, v = d->v;
+   long i, j, k, l, ct = 0, prec2;
+   gpmem_t av = avma;
+   FP_chk_fun chk = { &chk_reccoeff, &chk_reccoeff_init, NULL, 0 };
+   chk.data = (void*)d;
-   N   = degpol(nf[1]);
+   d->G = min(-10, -bit_accuracy(prec) >> 4);
-   G   = min(-10, -bit_accuracy(prec) >> 4);
+   eps = gpowgs(stoi(10), min(-8, (d->G >> 1)));
-   eps = gpowgs(stoi(10), min(-8, (G >> 1)));
    tB  = gpow(gmul2n(eps, N), gdivgs(gun, 1-N), DEFAULTPREC);
-   Bd    = gmin(B, tB);
+   Bd    = gceil(gmin(B, tB));
    p1    = gceil(gdiv(glog(Bd, DEFAULTPREC), dbltor(2.3026)));
    prec2 = max((prec << 1) - 2, (long)(itos(p1) * pariK1 + BIGDEFAULTPREC));
    nf2   = nfnewprec(nf, prec2);
-Line 2293  RecCoeff3(GEN nf, GEN beta, GEN B, long v, long prec)
+Line 2204  RecCoeff3(GEN nf, GEN beta, GEN B, long v, long prec)
 Line 2293  RecCoeff3(GEN nf, GEN beta, GEN B, long v, long prec)
 Line 2204  RecCoeff3(GEN nf, GEN beta, GEN B, long v, long prec)
    C2 = gdiv(B2, gsqr(eps));
    M = gmael(nf2, 5, 1);
+   d->M = M;
    A = cgetg(N+2, t_MAT);
    for (i = 1; i <= N+1; i++)
-Line 2318  RecCoeff3(GEN nf, GEN beta, GEN B, long v, long prec)
+Line 2230  RecCoeff3(GEN nf, GEN beta, GEN B, long v, long prec)
 Line 2318  RecCoeff3(GEN nf, GEN beta, GEN B, long v, long prec)
 Line 2230  RecCoeff3(GEN nf, GEN beta, GEN B, long v, long prec)
      }
    nB = mulsi(N+1, B2);
+   d->nB = nB;
+   cand = fincke_pohst(A, nB, 20000, 1, prec2, &chk);
-   data = new_chunk(8);
-   data[0] = (long)M;
-   data[1] = (long)beta;
-   data[2] = (long)B;
-   data[3] = v;
-   data[4] = G;
-   data[5] = N;
-   data[6] = (long)NULL;
-   data[7] = (long)nB;
-   chk = (FP_chk_fun*)new_chunk(sizeof(FP_chk_fun));
-   chk->f         = &chk_reccoeff;
-   chk->f_init    = &chk_reccoeff_init;
-   chk->f_post    = &chk_reccoeff_post;
-   chk->data      = data;
-   chk->skipfirst = 0;
-   cand = fincke_pohst(A, nB, 20000, 3, prec2, chk);
    if (!cand)
    {
      if (ct > 3) { avma = av; return NULL; }
-Line 2360  RecCoeff3(GEN nf, GEN beta, GEN B, long v, long prec)
+Line 2255  RecCoeff3(GEN nf, GEN beta, GEN B, long v, long prec)
 Line 2360  RecCoeff3(GEN nf, GEN beta, GEN B, long v, long prec)
 Line 2255  RecCoeff3(GEN nf, GEN beta, GEN B, long v, long prec)
    avma = av; return NULL;
  }
- /* Attempts to find a polynomial with coefficients in nf such that
+ /* Using linear dependance relations */
-    its coefficients are close to those of pol at the place v and
+ static GEN
+ RecCoeff2(GEN nf,  RC_data *d,  long prec)
+ {
+   long i, bacmin, bacmax;
+   gpmem_t av = avma, av2;
+   GEN vec, velt, p1, cand, M, plg, pol, cand2;
+   GEN beta = d->beta;
+   d->G = min(-20, -bit_accuracy(prec) >> 4);
+   M    = gmael(nf, 5, 1);
+   pol  = (GEN)nf[1];
+   vec  = gtrans((GEN)gtrans(M)[d->v]);
+   velt = (GEN)nf[7];
+   p1 = cgetg(2, t_VEC);
+   p1[1] = lneg(beta);
+   vec = concat(p1, vec);
+   p1[1] = zero;
+   velt = concat(p1, velt);
+   bacmin = (long)(.225 * bit_accuracy(prec));
+   bacmax = (long)(.315 * bit_accuracy(prec));
+   av2 = avma;
+   for (i = bacmax; i >= bacmin; i-=16)
+   {
+     p1 = lindep2(vec, i);
+     if (signe((GEN)p1[1]))
+     {
+       p1    = ground(gdiv(p1, (GEN)p1[1]));
+       cand  = gmodulcp(gmul(velt, gtrans(p1)), pol);
+       cand2 = algtobasis(nf, cand);
+       plg   = gmul(M, cand2);
+       if (TestOne(plg, d)) return gerepilecopy(av, cand);
+     }
+     avma = av2;
+   }
+   /* failure */
+   return RecCoeff3(nf,d,prec);
+ }
+ /* Attempts to find a polynomial with coefficients in nf such that
+    its coefficients are close to those of pol at the place v and
     less than B at all the other places */
  GEN
  RecCoeff(GEN nf,  GEN pol,  long v, long prec)
  {
-   long av = avma, j, md, G, cl = degpol(pol);
+   long j, md, G, cl = degpol(pol);
+   gpmem_t av = avma;
    GEN p1, beta;
+   RC_data d;
    /* if precision(pol) is too low, abort */
    for (j = 2; j <= cl+1; j++)
-Line 2380  RecCoeff(GEN nf,  GEN pol,  long v, long prec)
+Line 2323  RecCoeff(GEN nf,  GEN pol,  long v, long prec)
 Line 2380  RecCoeff(GEN nf,  GEN pol,  long v, long prec)
 Line 2323  RecCoeff(GEN nf,  GEN pol,  long v, long prec)
    md = cl/2;
    pol = dummycopy(pol);
+   d.N = degpol(nf[1]);
+   d.v = v;
    for (j = 1; j <= cl; j++)
-   {
+   { /* start with the coefficients in the middle,
-     /* start with the coefficients in the middle,
         since they are the harder to recognize! */
      long cf = md + (j%2? j/2: -j/2);
      GEN bound = binome(stoi(cl), cf);
      bound = shifti(bound, cl - cf);
-     if (DEBUGLEVEL > 1) fprintferr("In RecCoeff with cf = %ld and B = %Z\n",
+     if (DEBUGLEVEL > 1)
-                                    cf, bound);
+       fprintferr("In RecCoeff with cf = %ld and B = %Z\n", cf, bound);
      beta = greal((GEN)pol[cf+2]);
-     p1 = RecCoeff2(nf, beta, bound, v, prec);
+     d.beta = beta;
-     if (!p1)
+     d.B    = bound;
-     {
+     if (! (p1 = RecCoeff2(nf, &d, prec)) ) return NULL;
-       p1 = RecCoeff3(nf, beta, bound, v, prec);
-       if (!p1) return NULL;
-     }
      pol[cf+2] = (long)p1;
    }
    pol[cl+2] = un;
-Line 2406  RecCoeff(GEN nf,  GEN pol,  long v, long prec)
+Line 2348  RecCoeff(GEN nf,  GEN pol,  long v, long prec)
 Line 2406  RecCoeff(GEN nf,  GEN pol,  long v, long prec)
 Line 2348  RecCoeff(GEN nf,  GEN pol,  long v, long prec)
  }
  /*******************************************************************/
- /*******************************************************************/
  /*                                                                 */
- /*                   Computation of class fields of                */
+ /*     Class fields of real quadratic fields using Stark units     */
- /*                real quadratic fields using Stark units          */
  /*                                                                 */
  /*******************************************************************/
- /*******************************************************************/
- /* compute the coefficients an for the quadratic case */
+ /* an[q * i] *= chi for all (i,p)=1 */
- static int***
+ static void
- computean(GEN dtcr,  long nmax, long prec)
+ an_mul(int **an, long p, long q, long n, long deg, GEN chi, int **reduc)
  {
-   long i, j, cl, q, cp, al, v1, v2, v, fldiv, av, av1;
+   gpmem_t av;
-   int ***matan, ***reduc;
+   long c,i;
-   GEN bnf, ideal, dk, degs, idno, p1, prime, chi, qg, chi1, chi2;
+   int *T;
-   GEN chi11, chi12, bnr, pr, pr1, pr2, xray, xray1, xray2, dataray;
-   byteptr dp = diffptr;
+   if (gcmp1(chi)) return;
    av = avma;
+   T = (int*)new_chunk(deg); Polmod2Coeff(T,chi, deg);
+   for (c = 1, i = q; i <= n; i += q, c++)
+     if (c == p) c = 0; else MulCoeff(an[i], T, reduc, deg);
+   avma = av;
+ }
+ /* an[q * i] = 0 for all (i,p)=1 */
+ static void
+ an_set0_coprime(int **an, long p, long q, long n, long deg)
+ {
+   long c,i;
+   for (c = 1, i = q; i <= n; i += q, c++)
+     if (c == p) c = 0; else _0toCoeff(an[i], deg);
+ }
+ /* an[q * i] = 0 for all i */
+ static void
+ an_set0(int **an, long p, long n, long deg)
+ {
+   long i;
+   for (i = p; i <= n; i += p) _0toCoeff(an[i], deg);
+ }
-   cl = lg(dtcr) - 1;
+ /* compute the coefficients an for the quadratic case */
-   degs = GetDeg(dtcr);
+ static int**
+ computean(GEN dtcr, LISTray *R, long n, long deg)
+ {
+   gpmem_t av = avma, av2;
+   long i, p, q, condZ, l;
+   int **an, **reduc;
+   GEN L, CHI, chi, chi1, ray;
+   CHI_t C;
-   matan = InitMatAn(cl, nmax, degs, 1);
+   CHI = (GEN)dtcr[5]; init_CHI_alg(&C, CHI);
-   reduc = InitReduction(dtcr, degs);
+   condZ= R->condZ;
-   bnr = gmael(dtcr, 1, 4); bnf = (GEN)bnr[1];
+   an = InitMatAn(n, deg, 1);
-   dataray = InitGetRay(bnr, nmax);
+   reduc = InitReduction(CHI, deg);
+   av2 = avma;
-   ideal = gmael3(bnr, 2, 1, 1);
+   /* all pr | p divide cond */
-   idno  = idealnorm(bnf, ideal);
+   L = R->L0; l = lg(L);
-   dk = gmael(bnf, 7, 3);
+   for (i=1; i<l; i++) an_set0(an,L[i],n,deg);
-   prime = stoi(2);
+   /* 1 prime of degree 2 */
-   dp++;
+   L = R->L2; l = lg(L);
+   for (i=1; i<l; i++, avma = av2)
-   av1 = avma;
-   chi = chi1 = chi2 = NULL; /* gcc -Wall */
-   while (*dp && prime[2] <= nmax)
    {
-     qg = prime;
+     p = L[i];
-     al = 1;
+     if (condZ == 1) chi = C.val[0]; /* 1 */
+     else
+     {
+       ray = R->rayZ[p % condZ];
+       chi  = EvalChar(&C, ray);
+     }
+     chi1 = chi;
+     for (q=p;;)
+     {
+       an_set0_coprime(an, p,q,n,deg); /* v_p(q) odd */
+       if (! (q = next_pow(q,p, n)) ) break;
-     switch (krogs(dk, prime[2]))
+       an_mul(an,p,q,n,deg,chi,reduc);
+       if (! (q = next_pow(q,p, n)) ) break;
+       chi = gmul(chi, chi1);
+     }
+   }
+   /* 1 prime of degree 1 */
+   L = R->L1; l = lg(L);
+   for (i=1; i<l; i++, avma = av2)
+   {
+     p = L[i]; ray = R->L1ray[i];
+     chi  = EvalChar(&C, ray);
+     chi1 = chi;
+     for(q=p;;)
      {
-       /* prime is inert */
+       an_mul(an,p,q,n,deg,chi,reduc);
-       case -1:
+       if (! (q = next_pow(q,p, n)) ) break;
-         fldiv = divise(idno, prime);
+       chi = gmul(chi, chi1);
+     }
+   }
-         if (!fldiv)
+   /* 2 primes of degree 1 */
-         {
+   L = R->L11; l = lg(L);
-           xray = GetRay(bnr, dataray, prime, prec);
+   for (i=1; i<l; i++, avma = av2)
-           chi  = chiideal(dtcr, xray, 1);
+   {
-           chi1 = dummycopy(chi);
+     GEN ray1, ray2, chi11, chi12, chi2;
-         }
-         while(cmpis(qg, nmax) <= 0)
+     p = L[i]; ray1 = R->L11ray[i]; /* use pr1 pr2 = (p) */
-         {
+     if (condZ == 1)
-           q = qg[2];
+       ray2 = gneg(ray1);
+     else
+       ray2 = gsub(R->rayZ[p % condZ],  ray1);
+     chi11 = EvalChar(&C, ray1);
+     chi12 = EvalChar(&C, ray2);
-           for (cp = 1, i = q; i <= nmax; i += q, cp++)
+     chi1 = gadd(chi11, chi12);
-             if(cp % prime[2])
+     chi2 = chi12;
-             {
+     for(q=p;;)
-               if (fldiv || al%2)
+     {
-                 for (j = 1; j <= cl; j++)
+       an_mul(an,p,q,n,deg,chi1,reduc);
-                   _0toCoeff(matan[j][i], degs[j]);
+       if (! (q = next_pow(q,p, n)) ) break;
-               else
+       chi2 = gmul(chi2, chi12);
-                 for (j = 1; j <= cl; j++)
+       chi1 = gadd(chi2, gmul(chi1, chi11));
-                   MulPolmodCoeff((GEN)chi[j], matan[j][i], reduc[j], degs[j]);
-             }
-           qg = mulsi(q, prime);
-           al++;
-           if (al%2 && !fldiv)
-             for (j = 1; j <= cl; j++)
-               chi[j] = lmul((GEN)chi[j], (GEN)chi1[j]);
-         }
-         break;
-     /* prime is ramified */
-     case 0:
-       fldiv = divise(idno, prime);
-       if (!fldiv)
-       {
-         pr   = (GEN)primedec(bnf, prime)[1];
-         xray = GetRay(bnr, dataray, pr, prec);
-         chi  = chiideal(dtcr, xray, 1);
-         chi2 = dummycopy(chi);
-       }
-       while(cmpis(qg, nmax) <= 0)
-       {
-         q = qg[2];
-         for (cp = 1, i = q; i <= nmax; i += q, cp++)
-           if(cp % prime[2])
-           {
-             if (fldiv)
-               for(j = 1; j <= cl; j++)
-                 _0toCoeff(matan[j][i], degs[j]);
-             else
-             {
-               for(j = 1; j <= cl; j++)
-                 MulPolmodCoeff((GEN)chi[j], matan[j][i], reduc[j], degs[j]);
-             }
-           }
-         qg = mulsi(q, prime);
-         al++;
-         if (cmpis(qg, nmax) <= 0 && !fldiv)
-           for (j = 1; j <= cl; j++)
-             chi[j] = lmul((GEN)chi[j], (GEN)chi2[j]);
-       }
-       break;
-     /* prime is split */
-     default: /* case 1: */
-       p1  = primedec(bnf, prime);
-       pr1 = (GEN)p1[1];
-       pr2 = (GEN)p1[2];
-       v1 = idealval(bnf, ideal, pr1);
-       v2 = idealval(bnf, ideal, pr2);
-       if (v1 + v2 == 0)
-       {
-         xray1 = GetRay(bnr, dataray, pr1, prec);
-         xray2 = GetRay(bnr, dataray, pr2, prec);
-         chi11 = chiideal(dtcr, xray1, 1);
-         chi12 = chiideal(dtcr, xray2, 1);
-         chi1 = gadd(chi11, chi12);
-         chi2 = dummycopy(chi12);
-         while(cmpis(qg, nmax) <= 0)
-         {
-           q = qg[2];
-           for (cp = 1, i = q; i <= nmax; i += q, cp++)
-             if(cp % prime[2])
-               for(j = 1; j <= cl; j++)
-                 MulPolmodCoeff((GEN)chi1[j], matan[j][i], reduc[j], degs[j]);
-           qg = mulsi(q, prime);
-           al++;
-           if(cmpis(qg, nmax) <= 0)
-             for (j = 1; j <= cl; j++)
-             {
-               chi2[j] = lmul((GEN)chi2[j], (GEN)chi12[j]);
-               chi1[j] = ladd((GEN)chi2[j], gmul((GEN)chi1[j], (GEN)chi11[j]));
-             }
-         }
-       }
-       else
-       {
-         if (v1) { v  = v2; pr = pr2; } else { v  = v1; pr = pr1; }
-         if (v == 0)
-         {
-           xray = GetRay(bnr, dataray, pr, prec);
-           chi1 = chiideal(dtcr, xray, 1);
-           chi  = gcopy(chi1);
-         }
-         while(cmpis(qg, nmax) <= 0)
-         {
-           q = qg[2];
-           for (cp = 1, i = q; i <= nmax; i += q, cp++)
-             if(cp % prime[2])
-             {
-               if (v)
-                 for (j = 1; j <= cl; j++)
-                   _0toCoeff(matan[j][i], degs[j]);
-               else
-                 for (j = 1; j <= cl; j++)
-                   MulPolmodCoeff((GEN)chi[j], matan[j][i], reduc[j], degs[j]);
-             }
-           qg = mulii(qg, prime);
-           al++;
-           if (!v && (cmpis(qg, nmax) <= 0))
-             for (j = 1; j <= cl; j++)
-               chi[j] = lmul((GEN)chi[j], (GEN)chi1[j]);
-         }
-       }
-       break;
      }
-     prime[2] += (*dp++);
-     avma = av1;
    }
-   for (i = 1; i <= cl; i++)
+   CorrectCoeff(dtcr, an, reduc, n, deg);
-     CorrectCoeff((GEN)dtcr[i], matan[i], reduc[i], nmax, degs[i]);
+   FreeMat(reduc, deg-1);
+   avma = av; return an;
-   FreeMat(reduc);
-   avma = av; return matan;
  }
  /* compute S and T for the quadratic case */
  static GEN
- QuadGetST(GEN data, long prec)
+ QuadGetST(GEN dataCR, GEN vChar, long prec)
  {
-   long av = avma, n, j, nmax, cl, av1, av2, k;
+   const long cl  = lg(dataCR) - 1;
-   int ***matan;
+   gpmem_t av, av1, av2;
-   GEN nn, C, p1, p2, c2, cexp, cn, v, veclprime2, veclprime1;
+   long ncond, n, j, k, nmax;
-   GEN dtcr, cond, rep, an, cf, degs, veint1;
+   GEN rep, N0, C, T, S, cf, cfh, an, degs;
+   LISTray LIST;
-   dtcr     = (GEN)data[5];
+   /* allocate memory for answer */
-   cl       = lg(dtcr) - 1;
-   degs     = GetDeg(dtcr);
-   cf   = gmul2n(mpsqrt(mppi(prec)), 1);
-   C    = cgetg(cl+1, t_VEC);
-   cond = cgetg(cl+1, t_VEC);
-   c2   = cgetg(cl + 1, t_VEC);
-   nn   = new_chunk(cl+1);
-   nmax = 0;
-   for (j = 1; j <= cl; j++)
-   {
-     C[j]    = mael(dtcr, j, 2);
-     c2[j]   = ldivsg(2, (GEN)C[j]);
-     cond[j] = mael(dtcr, j, 7);
-     nn[j]   = (long)(bit_accuracy(prec) * gtodouble((GEN)C[j]) * 0.35);
-     nmax  = max(nmax, nn[j]);
-   }
-   if (nmax >= VERYBIGINT)
-     err(talker, "Too many coefficients (%ld) in QuadGetST: computation impossible", nmax);
-   if (DEBUGLEVEL >= 2)
-     fprintferr("nmax = %ld\n", nmax);
-   /* compute the coefficients */
-   matan = computean(dtcr, nmax, prec);
-   if (DEBUGLEVEL) msgtimer("Compute an");
-   /* allocate memory for the answer */
    rep = cgetg(3, t_VEC);
+   S = cgetg(cl+1, t_VEC); rep[1] = (long)S;
-   /* allocate memory for veclprime1 */
+   T = cgetg(cl+1, t_VEC); rep[2] = (long)T;
-   veclprime1 = cgetg(cl + 1, t_VEC);
    for (j = 1; j <= cl; j++)
    {
-     v = cgetg(3, t_COMPLEX);
+     S[j] = (long)cgetc(prec);
-     v[1] = lgetr(prec);
+     T[j] = (long)cgetc(prec);
-     v[2] = lgetr(prec); gaffect(gzero, v);
-     veclprime1[j] = (long)v;
    }
+   av = avma;
-   av1 = avma;
+   /* initializations */
-   cn = cgetr(prec);
+   degs = GetDeg(dataCR);
-   p1 = gmul2n(cf, -1);
+   ncond = lg(vChar)-1;
+   C    = cgetg(ncond+1, t_VEC);
-   /* compute veclprime1 */
+   N0   = cgetg(ncond+1, t_VECSMALL);
-   for (j = 1; j <= cl; j++)
+   nmax = 0;
+   for (j = 1; j <= ncond; j++)
    {
-     long n0 = 0;
+     C[j]  = mael(dataCR, mael(vChar,j,1), 2);
-     p2 = gmael3(dtcr, j, 5, 2);
+     N0[j] = (long)(bit_accuracy(prec) * gtodouble((GEN)C[j]) * 0.35);
-     cexp = gexp(gneg_i((GEN)c2[j]), prec);
+     if (nmax < N0[j]) nmax = N0[j];
-     av2 = avma; affsr(1, cn); v = (GEN)veclprime1[j];
-     for (n = 1; n <= nn[j]; n++)
-       if ( (an = EvalCoeff(p2, matan[j][n], degs[j])) )
-       {
-         affrr(gmul(cn, gpowgs(cexp, n - n0)), cn);
-         n0 = n;
-         gaffect(gadd(v, gmul(divrs(cn,n), an)), v);
-         avma = av2;
-       }
-     gaffect(gmul(p1, gmul(v, (GEN)C[j])), v);
-     avma = av2;
    }
-   avma = av1;
+   if ((ulong)nmax > maxprime())
-   rep[1] = (long)veclprime1;
+     err(talker, "Not enough precomputed primes (need all p <= %ld)", nmax);
-   if (DEBUGLEVEL) msgtimer("Compute V1");
+   if (DEBUGLEVEL>1) fprintferr("nmax = %ld\n", nmax);
+   InitPrimesQuad(gmael(dataCR,1,4), nmax, &LIST);
-   /* allocate memory for veclprime2 */
+   cf  = gmul2n(mpsqrt(mppi(prec)), 1);
-   veclprime2 = cgetg(cl + 1, t_VEC);
+   cfh = gmul2n(cf, -1);
-   for (j = 1; j <= cl; j++)
-   {
-     v = cgetg(3, t_COMPLEX);
-     v[1] = lgetr(prec);
-     v[2] = lgetr(prec); gaffect(gzero, v);
-     veclprime2[j] = (long)v;
-   }
-   /* compute f1(C/n) */
    av1 = avma;
+   /* loop over conductors */
-   veint1 = cgetg(cl + 1, t_VEC);
+   for (j = 1; j <= ncond; j++)
-   for (j = 1; j <= cl; j++)
    {
-     p1 = NULL;
+     const GEN c1 = (GEN)C[j], c2 = divsr(2,c1), cexp = mpexp(gneg(c2));
-     for (k = 1; k < j; k++)
+     const GEN LChar = (GEN)vChar[j];
-       if (gegal((GEN)cond[j], (GEN)cond[k])) { p1 = (GEN)veint1[k]; break; }
+     const long nChar = lg(LChar)-1, NN = N0[j];
-     if (p1 == NULL)
+     GEN veint1, vcn = cgetg(NN+1, t_VEC);
-     {
-       p1 = veceint1((GEN)c2[j], stoi(nn[j]), prec);
-       veint1[j] = (long)p1;
-     }
-     av2 = avma; p2 = gmael3(dtcr, j, 5, 2);
-     v = (GEN)veclprime2[j];
-     for (n = 1; n <= nn[j]; n++)
-       if ( (an = EvalCoeff(p2, matan[j][n], degs[j])) )
-       {
-         gaffect(gadd(v, gmul((GEN)p1[n], an)), v);
-         avma = av2;
-       }
-     gaffect(gmul(cf, gconj(v)), v);
-     avma = av2;
-   }
-   avma = av1;
-   rep[2] = (long)veclprime2;
-   if (DEBUGLEVEL) msgtimer("Compute V2");
-   FreeMat(matan); return gerepileupto(av, rep);
- }
- #if 0
+     if (DEBUGLEVEL>1)
- /* recover a quadratic integer by an exhaustive search */
+       fprintferr("* conductor no %ld/%ld (N = %ld)\n\tInit: ", j,ncond,NN);
- static GEN
+     veint1 = veceint1(c2, stoi(NN), prec);
- recbeta2(GEN nf,  GEN beta,  GEN bound,  long prec)
+     vcn[1] = (long)cexp;
- {
+     for (n=2; n<=NN; n++) vcn[n] = lmulrr((GEN)vcn[n-1], cexp);
-   long av = avma, av2, tetpil, i, range, G, e, m;
+     av2 = avma;
-   GEN om, om1, om2, dom, p1, a, b, rom, bom2, *gptr[2];
+     for (n=2; n<=NN; n++, avma = av2)
+       affrr(divrs((GEN)vcn[n],n), (GEN)vcn[n]);
-   G = min( - 20, - bit_accuracy(prec) >> 4);
+     for (k = 1; k <= nChar; k++)
-   if (DEBUGLEVEL > 3)
-     fprintferr("\n Precision needed: %ld", G);
-   om  = gmael(nf, 7, 2);
-   rom = (GEN)nf[6];
-   om1 = poleval(om, (GEN)rom[2]);
-   om2 = poleval(om, (GEN)rom[1]);
-   dom = subrr(om1, om2);
-   /* b will run from b to b + range */
-   p1 = gaddgs(gmul2n(gceil(absr(divir(bound, dom))), 1), 2);
-   range = VERYBIGINT;
-   if (cmpis(p1,  VERYBIGINT) < 0)
-     range = itos(p1);
-   av2 = avma;
-   b = gdiv(gsub(bound, beta), dom);
-   if (gsigne(b) < 0)
-     b = subis(negi(gcvtoi(gneg_i(b), &e)), 1);
-   else
-     b=gcvtoi(b, &e);
-   if (e > 0)  /* precision is lost in truncation */
-   {
-     avma = av;
-     return NULL;
-   }
-   bom2 = mulir(b, om2);
-   m = 0;
-   for (i = 0; i <= range; i++)
-   {
-     /* for every b,  we construct a and test it */
-     a = grndtoi(gsub(beta, bom2), &e);
-     if (e > 0) /* precision is lost in truncation */
      {
-       avma = av;
+       const long t = LChar[k], d = degs[t];
-       return NULL;
+       const GEN z = gmael3(dataCR, t, 5, 2);
-     }
+       GEN p1 = gzero, p2 = gzero;
+       int **matan;
+       long c = 0;
-     p1 = gsub(mpadd(a, bom2),  beta);
+       if (DEBUGLEVEL>1)
+         fprintferr("\tcharacter no: %ld (%ld/%ld)\n", t,k,nChar);
-     if ((DEBUGLEVEL > 3) && (expo(p1)<m))
+       matan = computean((GEN)dataCR[t], &LIST, NN, d);
-     {
+       for (n = 1; n <= NN; n++)
-       m = expo(p1);
+         if ((an = EvalCoeff(z, matan[n], d)))
-       fprintferr("\n Precision found: %ld", expo(p1));
+         {
+           p1 = gadd(p1, gmul(an, (GEN)vcn[n]));
+           p2 = gadd(p2, gmul(an, (GEN)veint1[n]));
+           if (++c == 256)
+           { GEN *gptr[2]; gptr[0]=&p1; gptr[1]=&p2;
+             gerepilemany(av2,gptr,2); c = 0;
+           }
+         }
+       gaffect(gmul(cfh, gmul(p1,c1)), (GEN)S[t]);
+       gaffect(gmul(cf,  gconj(p2)),   (GEN)T[t]);
+       FreeMat(matan,NN); avma = av2;
      }
+     if (DEBUGLEVEL>1) fprintferr("\n");
-     if (gcmp0(p1) || (expo(p1) < G))  /* result found */
+     avma = av1;
-     {
-       p1 = gadd(a, gmul(b, om));
-       return gerepileupto(av, gmodulcp(p1, (GEN)nf[1]));
-     }
-     tetpil = avma;
-     b    = gaddgs(b, 1);
-     bom2 = gadd(bom2, om2);
-     gptr[0] = &b;
-     gptr[1] = &bom2;
-     gerepilemanysp(av2, tetpil, gptr, 2);
    }
+   if (DEBUGLEVEL) msgtimer("S & T");
-   /* if it fails... */
+   avma = av; return rep;
-   return NULL;
  }
- #endif
+ typedef struct {
+   long cl;
+   GEN dkpow;
+ } DH_t;
  /* return 1 if the absolute polynomial pol (over Q) defines the
     Hilbert class field of the real quadratic field bnf */
- int
+ static GEN
- define_hilbert(GEN bnf, GEN pol)
+ define_hilbert(void *S, GEN pol)
  {
-   long cl;
+   DH_t *T = (DH_t*)S;
-   GEN dk;
+   GEN d = modulargcd(derivpol(pol), pol);
-   cl = itos(gmael3(bnf, 8, 1, 1));
+   if (degpol(pol) != T->cl + degpol(d)) return NULL;
-   dk = gmael(bnf, 7, 3);
+   pol = gdivexact(pol, d);
+   return (T->cl & 1 || !egalii(smalldiscf(pol), T->dkpow))? pol: NULL;
-   if (degpol(pol) == cl)
-     if ((cl%2) || !egalii(discf(pol), gpowgs(dk,cl>>1))) return 1;
-   return 0;
  }
  /* let polrel define Hk/k,  find L/Q such that Hk=Lk and L and k are
     disjoint */
  static GEN
- makescind(GEN bnf, GEN polabs, long cl, long prec)
+ makescind(GEN nf, GEN polrel, long cl)
  {
-   long av = avma, i, l;
+   long i, l;
-   GEN pol, p1, nf2, dabs, dk, bas;
+   gpmem_t av = avma;
+   GEN pol, polabs, L, BAS, nf2, dabs, dk, bas;
+   DH_t T;
+   FP_chk_fun CHECK;
-   /* check the result (a little): signature and discriminant */
+   BAS = rnfpolredabs(nf, polrel, nf_ABSOLUTE|nf_ADDZK);
-   bas = allbase4(polabs,0,&dabs,NULL);
+   polabs = (GEN)BAS[1];
-   dk  = gmael(bnf,7,3);
+   bas    = (GEN)BAS[2];
+   dabs = gmul(ZX_disc(polabs), gsqr(det2(vecpol_to_mat(bas, 2*cl))));
+   /* check result (a little): signature and discriminant */
+   dk  = (GEN)nf[3];
    if (!egalii(dabs, gpowgs(dk,cl)) || sturm(polabs) != 2*cl)
      err(bugparier, "quadhilbert");
    /* attempt to find the subfields using polred */
-   p1 = cgetg(3,t_VEC); p1[1]=(long)polabs; p1[2]=(long)bas;
+   T.cl = cl;
-   pol = polredfirstpol(p1, (prec<<1) - 2, &define_hilbert, bnf);
+   T.dkpow = (cl & 1) ? NULL: gpowgs(dk, cl>>1);
+   CHECK.f = &define_hilbert;
+   CHECK.data = (void*)&T;
+   pol = polredfirstpol(BAS, 0, &CHECK);
    if (DEBUGLEVEL) msgtimer("polred");
    if (!pol)
    {
-     nf2 = nfinit0(polabs, 1, prec);
+     nf2 = nfinit0(BAS, 0, DEFAULTPREC);
-     p1  = subfields(nf2, stoi(cl));
+     L  = subfields(nf2, stoi(cl));
-     l = lg(p1);
+     l = lg(L);
      if (DEBUGLEVEL) msgtimer("subfields");
      for (i = 1; i < l; i++)
      {
-       pol = gmael(p1, i, 1);
+       pol = gmael(L, i, 1);
-       if ((cl%2) || !gegal(sqri(discf(pol)), (GEN)nf2[3])) break;
+       if (cl & 1 || !egalii(smalldiscf(pol), T.dkpow)) break;
      }
      if (i == l)
        for (i = 1; i < l; i++)
        {
-         pol = gmael(p1, i, 1);
+         pol = gmael(L, i, 1);
-         if (degpol(gcoeff(nffactor(bnf, pol), 1, 1)) == cl) break;
+         if (degpol(gcoeff(nffactor(nf, pol), 1, 1)) == cl) break;
        }
      if (i == l)
        err(bugparier, "makescind (no polynomial found)");
    }
-   pol = polredabs(pol, prec);
+   pol = polredabs0(pol, nf_PARTIALFACT);
    return gerepileupto(av, pol);
  }
-Line 2878  makescind(GEN bnf, GEN polabs, long cl, long prec)
+Line 2641  makescind(GEN bnf, GEN polabs, long cl, long prec)
 Line 2878  makescind(GEN bnf, GEN polabs, long cl, long prec)
 Line 2641  makescind(GEN bnf, GEN polabs, long cl, long prec)
  static GEN
  GenusField(GEN bnf, long prec)
  {
-   long hk, c, l, av = avma;
+   long hk, c, l;
-   GEN disc, quat, x2, pol, div, d;
+   gpmem_t av = avma;
+   GEN disc, x2, pol, div, d;
    hk   = itos(gmael3(bnf, 8, 1, 1));
    disc = gmael(bnf, 7, 3);
-   quat = stoi(4);
    x2   = gsqr(polx[0]);
-   if (gcmp0(modii(disc, quat))) disc = divii(disc, quat);
+   if (mod4(disc) == 0) disc = divis(disc, 4);
    div = divisors(disc);
-   c = 1;
    l = 0;
-   pol = NULL; /* gcc -Wall */
+   pol = NULL;
-   while(l < hk)
+   for (c = 2; l < hk; c++) /* skip c = 1 : div[1] = gun */
    {
-     c++;
      d = (GEN)div[c];
+     if (mod4(d) == 1)
-     if (gcmp1(modii(d, quat)))
      {
-       if (!l)
+       GEN t = gsub(x2, d);
-         pol = gsub(x2, d);
+       if (!pol)
+         pol = t;
        else
-         pol=(GEN)compositum(pol, gsub(x2, d))[1];
+         pol = (GEN)compositum(pol, t)[1];
        l = degpol(pol);
      }
    }
-   return gerepileupto(av, polredabs(pol, prec));
+   return gerepileupto(av, polredabs0(pol, nf_PARTIALFACT));
  }
  /* if flag = 0 returns the reduced polynomial,  flag = 1 returns the
-Line 2919  GenusField(GEN bnf, long prec)
+Line 2679  GenusField(GEN bnf, long prec)
 Line 2919  GenusField(GEN bnf, long prec)
 Line 2679  GenusField(GEN bnf, long prec)
  static GEN
  AllStark(GEN data,  GEN nf,  long flag,  long newprec)
  {
-   long cl, i, j, cpt = 0, av, av2, N, h, v, n, bnd = 300, sq = 1, r1, r2;
+   long cl, i, j, cpt = 0, N, h, v, n, bnd = 300, sq = 1, r1, r2;
-   int ***matan;
+   gpmem_t av, av2;
-   GEN p0, p1, p2, S, T, polrelnum, polrel, Lp, W, A, veczeta, sig, valchi;
+   int **matan;
-   GEN degs, ro, C, Cmax, dataCR, cond1, L1, *gptr[2], an, Pi;
+   GEN p0, p1, p2, S, T, polrelnum, polrel, Lp, W, A, veczeta, sig;
+   GEN vChar, degs, ro, C, Cmax, dataCR, cond1, L1, *gptr[2], an, Pi;
+   LISTray LIST;
    N     = degpol(nf[1]);
    r1    = nf_get_r1(nf);
    r2    = (N - r1)>>1;
    cond1 = gmael4(data, 1, 2, 1, 2);
    Pi    = mppi(newprec);
+   dataCR = (GEN)data[5];
    v = 1;
    while(gcmp1((GEN)cond1[v])) v++;
- LABDOUB:
-   av = avma;
-   dataCR = (GEN)data[5];
    cl = lg(dataCR)-1;
    degs = GetDeg(dataCR);
    h  = itos(gmul2n(det((GEN)data[2]), -1));
+ LABDOUB:
+   av = avma;
+   vChar = sortChars(dataCR, N==2);
+   W = ComputeAllArtinNumbers(dataCR, vChar, (flag >= 0), newprec);
    if (flag >= 0)
    {
-     /* compute S,T differently if nf is quadratic */
+     p1 = GetST(dataCR, vChar, newprec);
-     if (N == 2)
-       p1 = QuadGetST(data, newprec);
-     else
-       p1 = GetST(dataCR, newprec);
      S = (GEN)p1[1];
      T = (GEN)p1[2];
      Lp = cgetg(cl + 1, t_VEC);
      for (i = 1; i <= cl; i++)
-       Lp[i] = GetValue((GEN)dataCR[i], (GEN)S[i], (GEN)T[i], 0, 1, newprec)[2];
+       Lp[i] = GetValue((GEN)dataCR[i], (GEN)W[i], (GEN)S[i], (GEN)T[i],
+, 1, newprec)[2];
      if (DEBUGLEVEL) msgtimer("Compute W");
    }
    else
-   {
+   { /* compute a crude approximation of the result */
-     /* compute a crude approximation of the result */
      C = cgetg(cl + 1, t_VEC);
      for (i = 1; i <= cl; i++) C[i] = mael(dataCR, i, 2);
      Cmax = vecmax(C);
-     n = GetBoundN0(Cmax, r1, r2, newprec, 0);
+     n = GetBoundN0(Cmax, r1, r2, newprec);
      if (n > bnd) n = bnd;
      if (DEBUGLEVEL) fprintferr("nmax in QuickPol: %ld \n", n);
+     InitPrimes(gmael(dataCR,1,4), n, &LIST);
-     matan = ComputeCoeff(dataCR, n, newprec);
      p0 = cgetg(3, t_COMPLEX);
-     p0[1] = lgetr(newprec); affsr(0, (GEN)p0[1]);
+     p0[1] = (long)realzero(newprec);
-     p0[2] = lgetr(newprec); affsr(0, (GEN)p0[2]);
+     p0[2] = (long)realzero(newprec);
      L1 = cgetg(cl+1, t_VEC);
      /* we use the formulae L(1) = sum (an / n) */
      for (i = 1; i <= cl; i++)
      {
+       matan = ComputeCoeff((GEN)dataCR[i], &LIST, n, degs[i]);
        av2 = avma;
        p1 = p0; p2 = gmael3(dataCR, i, 5, 2);
        for (j = 1; j <= n; j++)
-         if ( (an = EvalCoeff(p2, matan[i][j], degs[i])) )
+         if ( (an = EvalCoeff(p2, matan[j], degs[i])) )
            p1 = gadd(p1, gdivgs(an, j));
        L1[i] = lpileupto(av2, p1);
+       FreeMat(matan, n);
      }
-     FreeMat(matan);
      p1 = gmul2n(gpowgs(mpsqrt(Pi), N - 2), 1);
      Lp = cgetg(cl+1, t_VEC);
      for (i = 1; i <= cl; i++)
      {
-       W = ComputeArtinNumber((GEN)dataCR[i], 1, newprec);
+       GEN dtcr = (GEN)dataCR[i], WW;
-       A = (GEN)ComputeAChi((GEN)dataCR[i], 0, newprec)[2];
+       A = (GEN)ComputeAChi(dtcr, 0, newprec)[2];
-       W = gmul((GEN)C[i], gmul(A, W));
+       WW= gmul((GEN)C[i], gmul(A, (GEN)W[i]));
-       Lp[i] = ldiv(gmul(W, gconj((GEN)L1[i])), p1);
+       Lp[i] = ldiv(gmul(WW, gconj((GEN)L1[i])), p1);
      }
    }
-Line 3010  LABDOUB:
+Line 2767  LABDOUB:
 Line 3010  LABDOUB:
 Line 2767  LABDOUB:
      GEN z = gzero;
      sig = (GEN)p1[i];
-     valchi = chiideal(dataCR, sig, 0);
      for (j = 1; j <= cl; j++)
      {
-       GEN p2 = greal(gmul((GEN)Lp[j], (GEN)valchi[j]));
+       GEN CHI = gmael(dataCR,j,5);
-       if (!gegal(gdeux, gmael3(dataCR, j, 5, 3)))
+       GEN val = ComputeImagebyChar(CHI, sig);
-         p2 = gmul2n(p2, 1); /* character not real */
+       GEN p2 = greal(gmul((GEN)Lp[j], val));
-       z = gadd(z,p2);
+       if (itos((GEN)CHI[3]) != 2) p2 = gmul2n(p2, 1); /* character not real */
+       z = gadd(z, p2);
      }
      veczeta[i] = ldivgs(z, 2 * h);
    }
    if (DEBUGLEVEL >= 2) fprintferr("zetavalues = %Z\n", veczeta);
-   if ((flag >=0) && (flag <= 3)) sq = 0;
+   if (flag >=0 && flag <= 3) sq = 0;
    ro = cgetg(h+1, t_VEC); /* roots */
    for (;;)
    {
-     if (!sq && (DEBUGLEVEL > 1))
+     if (DEBUGLEVEL > 1 && !sq)
        fprintferr("Checking the square-root of the Stark unit...\n");
      for (j = 1; j <= h; j++)
-Line 3043  LABDOUB:
+Line 2799  LABDOUB:
 Line 3043  LABDOUB:
 Line 2799  LABDOUB:
        if (DEBUGLEVEL >= 2) fprintferr("polrelnum = %Z\n", polrelnum);
        msgtimer("Compute %s", (flag < 0)? "quickpol": "polrelnum");
      }
      if (flag < 0)
        return gerepilecopy(av, polrelnum);
      /* we try to recognize this polynomial */
      polrel = RecCoeff(nf, polrelnum, v, newprec);
-Line 3066  LABDOUB:
+Line 2822  LABDOUB:
 Line 3066  LABDOUB:
 Line 2822  LABDOUB:
      if (DEBUGLEVEL) err(warnprec, "AllStark", newprec);
      nf = nfnewprec(nf, newprec);
-     data[5] = (long)CharNewPrec((GEN)data[5], nf, newprec);
+     dataCR = CharNewPrec(dataCR, nf, newprec);
-     gptr[0] = &data;
+     gptr[0] = &dataCR;
-     gptr[1] = &nf;
+     gptr[1] = &nf; gerepilemany(av, gptr, 2);
-     gerepilemany(av, gptr, 2);
      goto LABDOUB;
    }
-Line 3084  LABDOUB:
+Line 2839  LABDOUB:
 Line 3084  LABDOUB:
 Line 2839  LABDOUB:
      polrel = gmul(polrel, gpowgs(polx[n], h));
      polrel = gsubst(polrel, n, polx[0]);
      polrel = gmul(polrel, leading_term(polrel));
-     delete_var();
+     (void)delete_var();
    }
    if (DEBUGLEVEL >= 2) fprintferr("polrel = %Z\n", polrel);
    if (DEBUGLEVEL) msgtimer("Recpolnum");
    /* we want a reduced relative polynomial */
-   if (!flag) return gerepileupto(av, rnfpolredabs(nf, polrel, 0, newprec));
+   if (!flag) return gerepileupto(av, rnfpolredabs(nf, polrel, 0));
    /* we just want the polynomial computed */
    if (flag!=2) return gerepilecopy(av, polrel);
    /* we want a reduced absolute polynomial */
-   return gerepileupto(av, rnfpolredabs(nf, polrel, 2, newprec));
+   return gerepileupto(av, rnfpolredabs(nf, polrel, nf_ABSOLUTE));
  }
  /********************************************************************/
-Line 3109  LABDOUB:
+Line 2864  LABDOUB:
 Line 3109  LABDOUB:
 Line 2864  LABDOUB:
  GEN
  quadhilbertreal(GEN D, GEN flag, long prec)
  {
-   VOLATILE long av = avma, cl;
+   gpmem_t av = avma;
-   long newprec;
+   long cl, newprec;
-   VOLATILE GEN pol, bnf, bnr, dataC, bnrh, nf, exp;
+   GEN pol, bnf, bnr, dataC, bnrh, nf, exp;
-   void *catcherr = NULL;
-   if (DEBUGLEVEL) timer2();
+   (void)&prec; /* prevent longjmp clobbering it */
+   if (DEBUGLEVEL) (void)timer2();
    disable_dbg(0);
    /* quick computation of the class number */
-Line 3130  quadhilbertreal(GEN D, GEN flag, long prec)
+Line 2885  quadhilbertreal(GEN D, GEN flag, long prec)
 Line 3130  quadhilbertreal(GEN D, GEN flag, long prec)
 Line 2885  quadhilbertreal(GEN D, GEN flag, long prec)
    pol = quadpoly(D);
    setvarn(pol, fetch_var());
-  START:
+ START:
    /* compute the class group */
    bnf = bnfinit0(pol, 1, NULL, prec);
    nf  = (GEN)bnf[7];
    disable_dbg(-1);
    if (DEBUGLEVEL) msgtimer("Compute Cl(k)");
    /* if the exponent of the class group is 2, use rather Genus Field Theory */
    exp = gmael4(bnf, 8, 1, 2, 1);
-   if (gegal(exp, gdeux)) { delete_var(); return GenusField(bnf, prec); }
+   if (gegal(exp, gdeux)) { (void)delete_var(); return GenusField(bnf, prec); }
-   { /* catch precision problems (precision too small) */
+   CATCH(precer) {
-     jmp_buf env;
+     prec += EXTRA_PREC; pol = NULL;
-     if (setjmp(env))
+     err (warnprec, "quadhilbertreal", prec);
+   } TRY {
+     /* find the modulus defining N */
+     bnr   = buchrayinitgen(bnf, gun);
+     dataC = InitQuotient(bnr, gzero);
+     bnrh  = FindModulus(dataC, 1, &newprec, prec, gcmp0(flag)? 0: -10);
+     if (DEBUGLEVEL) msgtimer("FindModulus");
+     if (newprec > prec)
      {
-       prec += EXTRA_PREC;
+       if (DEBUGLEVEL >= 2) fprintferr("new precision: %ld\n", newprec);
-       err (warnprec, "quadhilbertreal", prec);
+       nf = nfnewprec(nf, newprec);
-       goto START;
      }
-     catcherr = err_catch(precer, env, NULL);
+     pol = AllStark(bnrh, nf, 1, newprec);
-   }
+   } ENDCATCH;
+   if (!pol) goto START;
-   /* find the modulus defining N */
+   pol = makescind(nf, pol, cl);
-   bnr   = buchrayinitgen(bnf, gun);
+   (void)delete_var();
-   dataC = InitQuotient(bnr, gzero);
+   return gerepileupto(av, pol);
-   bnrh  = FindModulus(dataC, 1, &newprec, prec, gcmp0(flag)? 0: -10);
+ }
-   if (DEBUGLEVEL) msgtimer("FindModulus");
+ static GEN
+ get_subgroup(GEN subgp, GEN cyc)
-   if (newprec > prec)
+ {
-   {
+   if (!subgp || gcmp0(subgp)) return cyc;
-     if (DEBUGLEVEL >= 2) fprintferr("new precision: %ld\n", newprec);
+   return gcmp1(denom(gauss(subgp, cyc)))? subgp: NULL;
-     nf = nfnewprec(nf, newprec);
-   }
-   /* use the generic function AllStark */
-   pol = AllStark(bnrh, nf, 2, newprec);
-   delete_var();
-   return gerepileupto(av, makescind(bnf, pol, cl, prec));
  }
  GEN
- bnrstark(GEN bnr,  GEN subgroup,  long flag,  long prec)
+ bnrstark(GEN bnr,  GEN subgrp,  long flag,  long prec)
  {
-   long cl, N, newprec, av = avma, bnd = 0;
+   long cl, N, newprec, bnd = 0;
+   gpmem_t av = avma;
    GEN bnf, dataS, p1, Mcyc, nf, data;
    if (flag >= 4)
    {
      bnd = -10;
      flag -= 4;
-Line 3188  bnrstark(GEN bnr,  GEN subgroup,  long flag,  long pre
+Line 2945  bnrstark(GEN bnr,  GEN subgroup,  long flag,  long pre
 Line 3188  bnrstark(GEN bnr,  GEN subgroup,  long flag,  long pre
 Line 2945  bnrstark(GEN bnr,  GEN subgroup,  long flag,  long pre
    /* check the bnr */
    checkbnrgen(bnr);
+   bnf  = checkbnf(bnr);
-   bnf  = (GEN)bnr[1];
+   nf   = checknf(bnf);
-   nf   = (GEN)bnf[7];
-   Mcyc = diagonal(gmael(bnr, 5, 2));
    N    = degpol(nf[1]);
-   if (N == 1)
+   if (N == 1) return galoissubcyclo(bnr, subgrp, 0, 0);
-     err(talker, "the ground field must be distinct from Q");
+   Mcyc = diagonal(gmael(bnr, 5, 2));
    /* check the bnf */
-   if (!varn(gmael(bnf, 7, 1)))
+   if (!varn(nf[1]))
      err(talker, "main variable in bnrstark must not be x");
-   if (cmpis(gmael3(bnf, 7, 2, 1), N))
+   if (nf_get_r2(nf))
      err(talker, "not a totally real ground base field in bnrstark");
-   /* check the subgroup */
+   /* check the subgrp */
-   if (gcmp0(subgroup))
+   if (! (subgrp = get_subgroup(subgrp,Mcyc)) )
-     subgroup = Mcyc;
+     err(talker, "incorrect subgrp in bnrstark");
-   else
-   {
-     p1 = gauss(subgroup, Mcyc);
-     if (!gcmp1(denom(p1)))
-       err(talker, "incorrect subgroup in bnrstark");
-   }
    /* compute bnr(conductor) */
-   p1       = conductor(bnr, subgroup, 2);
+   p1       = conductor(bnr, subgrp, 2);
    bnr      = (GEN)p1[2];
-   subgroup = (GEN)p1[3];
+   subgrp = (GEN)p1[3];
    /* check the class field */
    if (!gcmp0(gmael3(bnr, 2, 1, 2)))
      err(talker, "not a totally real class field in bnrstark");
-   cl = itos(det(subgroup));
+   cl = itos(det(subgrp));
    if (cl == 1) return polx[0];
-   timer2();
+   if (DEBUGLEVEL) (void)timer2();
    /* find a suitable extension N */
-   dataS = InitQuotient(bnr, subgroup);
+   dataS = InitQuotient(bnr, subgrp);
    data  = FindModulus(dataS, 1, &newprec, prec, bnd);
    if (newprec > prec)
    {
      if (DEBUGLEVEL >= 2) fprintferr("new precision: %ld\n", newprec);
      nf = nfnewprec(nf, newprec);
    }
    return gerepileupto(av, AllStark(data, nf, flag, newprec));
  }
- /* For each character of Cl(bnr)/sbgrp, compute L(1, chi) (or equivalently
+ /* For each character of Cl(bnr)/subgp, compute L(1, chi) (or equivalently
     the first non-zero term c(chi) of the expansion at s = 0). The binary
     digits of flag mean 1: if 0 then compute the term c(chi) and return
     [r(chi), c(chi)] where r(chi) is the order of L(s, chi) at s = 0,
-Line 3251  bnrstark(GEN bnr,  GEN subgroup,  long flag,  long pre
+Line 3001  bnrstark(GEN bnr,  GEN subgroup,  long flag,  long pre
 Line 3251  bnrstark(GEN bnr,  GEN subgroup,  long flag,  long pre
 Line 3001  bnrstark(GEN bnr,  GEN subgroup,  long flag,  long pre
     the modulus of bnr (and the infinite places), 3: return also the
     character */
  GEN
- bnrL1(GEN bnr, GEN sbgrp, long flag, long prec)
+ bnrL1(GEN bnr, GEN subgp, long flag, long prec)
  {
    GEN bnf, nf, cyc, Mcyc, p1, L1, chi, lchi, clchi, allCR, listCR, dataCR;
-   GEN S, T, rep, indCR, invCR, Qt;
+   GEN S, T, rep, indCR, invCR, Qt, vChar;
-   long N, cl, i, j, k, nc, lq, a, av = avma, ncc;
+   long N, cl, i, j, k, nc, lq, a, ncc;
+   gpmem_t av = avma;
-   checkbnr(bnr);
+   checkbnrgen(bnr);
    bnf  = (GEN)bnr[1];
    nf   = (GEN)bnf[7];
    cyc  = gmael(bnr, 5, 2);
-Line 3265  bnrL1(GEN bnr, GEN sbgrp, long flag, long prec)
+Line 3016  bnrL1(GEN bnr, GEN sbgrp, long flag, long prec)
 Line 3265  bnrL1(GEN bnr, GEN sbgrp, long flag, long prec)
 Line 3016  bnrL1(GEN bnr, GEN sbgrp, long flag, long prec)
    ncc  = lg(cyc) - 1;
    N    = degpol(nf[1]);
-   if (N == 1)
+   if (N == 1) err(talker, "the ground field must be distinct from Q");
-     err(talker, "the ground field must be distinct from Q");
+   if (flag < 0 || flag > 8) err(flagerr,"bnrL1");
-   if ((flag < 0) || (flag > 8))
-     err(flagerr,"bnrL1");
-   /* check the bnr */
-   checkbnrgen(bnr);
    /* compute bnr(conductor) */
    if (!(flag & 2))
    {
-Line 3284  bnrL1(GEN bnr, GEN sbgrp, long flag, long prec)
+Line 3029  bnrL1(GEN bnr, GEN sbgrp, long flag, long prec)
 Line 3284  bnrL1(GEN bnr, GEN sbgrp, long flag, long prec)
 Line 3029  bnrL1(GEN bnr, GEN sbgrp, long flag, long prec)
    }
    /* check the subgroup */
-   if (gcmp0(sbgrp))
+   if (! (subgp = get_subgroup(subgp,Mcyc)) )
-     sbgrp = Mcyc;
+     err(talker, "incorrect subgroup in bnrL1");
-   else
-   {
-     if (lg(sbgrp) != ncc+1)
-       err(talker, "incorrect subgroup in bnrL1");
-     p1 = gauss(sbgrp, Mcyc);
-     if (!gcmp1(denom(p1)))
-       err(talker, "incorrect subgroup in bnrL1");
-   }
-   cl = labs(itos(det(sbgrp)));
+   cl = labs(itos(det(subgp)));
-   Qt = InitQuotient0(Mcyc, sbgrp);
+   Qt = InitQuotient0(Mcyc, subgp);
    lq = lg((GEN)Qt[2]) - 1;
-   /* compute all the characters */
+   /* compute all characters */
-   allCR = FindEltofGroup(cl, (GEN)Qt[2]);
+   allCR = EltsOfGroup(cl, (GEN)Qt[2]);
    /* make a list of all non-trivial characters modulo conjugation */
    listCR = cgetg(cl, t_VEC);
-Line 3316  bnrL1(GEN bnr, GEN sbgrp, long flag, long prec)
+Line 3052  bnrL1(GEN bnr, GEN sbgrp, long flag, long prec)
 Line 3316  bnrL1(GEN bnr, GEN sbgrp, long flag, long prec)
 Line 3052  bnrL1(GEN bnr, GEN sbgrp, long flag, long prec)
      /* lift the character to a character on Cl(bnr) */
      lchi = cgetg(ncc + 1, t_VEC);
      for (j = 1; j <= ncc; j++)
      {
        p1 = gzero;
        for (k = 1; k <= lq; k++)
          p1 = gadd(p1, gdiv(mulii(gmael3(Qt, 3, j, k), (GEN)chi[k]),
                             gmael(Qt, 2, k)));
        lchi[j] = lmodii(gmul(p1, (GEN)cyc[j]), (GEN)cyc[j]);
      }
-Line 3355  bnrL1(GEN bnr, GEN sbgrp, long flag, long prec)
+Line 3091  bnrL1(GEN bnr, GEN sbgrp, long flag, long prec)
 Line 3355  bnrL1(GEN bnr, GEN sbgrp, long flag, long prec)
 Line 3091  bnrL1(GEN bnr, GEN sbgrp, long flag, long prec)
    /* compute the data for these characters */
    dataCR = InitChar(bnr, listCR, prec);
-   p1 = GetST(dataCR, prec);
+   vChar= sortChars(dataCR,0);
+   p1 = GetST(dataCR, vChar, prec);
    S = (GEN)p1[1];
    T = (GEN)p1[2];
-Line 3369  bnrL1(GEN bnr, GEN sbgrp, long flag, long prec)
+Line 3105  bnrL1(GEN bnr, GEN sbgrp, long flag, long prec)
 Line 3369  bnrL1(GEN bnr, GEN sbgrp, long flag, long prec)
 Line 3105  bnrL1(GEN bnr, GEN sbgrp, long flag, long prec)
    {
      a = indCR[i];
      if (a > 0)
-       L1[i] = (long)GetValue((GEN)dataCR[a], (GEN)S[a], (GEN)T[a], flag & 1,
+       L1[i] = (long)GetValue((GEN)dataCR[a], NULL, (GEN)S[a], (GEN)T[a],
-                              flag & 2, prec);
+                              flag & 1, flag & 2, prec);
    }
    for (i = 1; i < cl; i++)

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>