Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - modules - stark.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.10.0 lcov report (development 19837-cc815bb) Lines: 1813 1949 93.0 %
Date: 2016-12-10 05:49:10 Functions: 126 127 99.2 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : /*******************************************************************/
      15             : /*                                                                 */
      16             : /*        COMPUTATION OF STARK UNITS OF TOTALLY REAL FIELDS        */
      17             : /*                                                                 */
      18             : /*******************************************************************/
      19             : #include "pari.h"
      20             : #include "paripriv.h"
      21             : 
      22             : static const long EXTRA_PREC = DEFAULTPREC-2;
      23             : 
      24             : /* ComputeCoeff */
      25             : typedef struct {
      26             :   GEN L0, L1, L11, L2; /* VECSMALL of p */
      27             :   GEN L1ray, L11ray; /* precomputed isprincipalray(pr), pr | p */
      28             :   GEN rayZ; /* precomputed isprincipalray(i), i < condZ */
      29             :   long condZ; /* generates cond(bnr) \cap Z, assumed small */
      30             : } LISTray;
      31             : 
      32             : /* Char evaluation */
      33             : typedef struct {
      34             :   long ord;
      35             :   GEN *val, chi;
      36             : } CHI_t;
      37             : 
      38             : /* RecCoeff */
      39             : typedef struct {
      40             :   GEN M, beta, B, U, nB;
      41             :   long v, G, N;
      42             : } RC_data;
      43             : 
      44             : /********************************************************************/
      45             : /*                    Miscellaneous functions                       */
      46             : /********************************************************************/
      47             : static GEN
      48       18963 : chi_get_c(GEN chi) { return gmael(chi,1,2); }
      49             : static GEN
      50       58338 : chi_get_gdeg(GEN chi) { return gmael(chi,1,1); }
      51             : static long
      52       58338 : chi_get_deg(GEN chi) { return itou(chi_get_gdeg(chi)); }
      53             : 
      54             : /* Compute the image of logelt by character chi, as a complex number */
      55             : static ulong
      56       14868 : CharEval_n(GEN chi, GEN logelt)
      57             : {
      58       14868 :   GEN gn = ZV_dotproduct(chi_get_c(chi), logelt);
      59       14868 :   return umodiu(gn, chi_get_deg(chi));
      60             : }
      61             : /* Compute the image of logelt by character chi, as a complex number */
      62             : static GEN
      63       14735 : CharEval(GEN chi, GEN logelt)
      64             : {
      65       14735 :   ulong n = CharEval_n(chi, logelt), d = chi_get_deg(chi);
      66       14735 :   long nn = Fl_center(n,d,d>>1);
      67       14735 :   GEN x = gel(chi,2);
      68       14735 :   x = gpowgs(x, labs(nn));
      69       14735 :   if (nn < 0) x = gconj(x);
      70       14735 :   return x;
      71             : }
      72             : 
      73             : /* return n such that C(elt) = z^n */
      74             : static ulong
      75      679394 : CHI_eval_n(CHI_t *C, GEN logelt)
      76             : {
      77      679394 :   GEN n = ZV_dotproduct(C->chi, logelt);
      78      679394 :   return umodiu(n, C->ord);
      79             : }
      80             : /* return C(elt) */
      81             : static GEN
      82      677602 : CHI_eval(CHI_t *C, GEN logelt)
      83             : {
      84      677602 :   return C->val[CHI_eval_n(C, logelt)];
      85             : }
      86             : 
      87             : static void
      88        4095 : init_CHI(CHI_t *c, GEN CHI, GEN z)
      89             : {
      90        4095 :   long i, d = chi_get_deg(CHI);
      91        4095 :   GEN *v = (GEN*)new_chunk(d);
      92        4095 :   v[0] = gen_1;
      93        4095 :   if (d != 1)
      94             :   {
      95        4095 :     v[1] = z;
      96        4095 :     for (i=2; i<d; i++) v[i] = gmul(v[i-1], z);
      97             :   }
      98        4095 :   c->chi = chi_get_c(CHI);
      99        4095 :   c->ord = d;
     100        4095 :   c->val = v;
     101        4095 : }
     102             : /* as t_POLMOD */
     103             : static void
     104        2562 : init_CHI_alg(CHI_t *c, GEN CHI) {
     105        2562 :   long d = chi_get_deg(CHI);
     106             :   GEN z;
     107        2562 :   switch(d)
     108             :   {
     109           0 :     case 1: z = gen_1; break;
     110        1036 :     case 2: z = gen_m1; break;
     111        1526 :     default: z = mkpolmod(pol_x(0), polcyclo(d,0));
     112             :   }
     113        2562 :   init_CHI(c,CHI, z);
     114        2562 : }
     115             : /* as t_COMPLEX */
     116             : static void
     117        1533 : init_CHI_C(CHI_t *c, GEN CHI) {
     118        1533 :   init_CHI(c,CHI, gel(CHI,2));
     119        1533 : }
     120             : 
     121             : typedef struct {
     122             :   long r; /* rank = lg(gen) */
     123             :   GEN j; /* current elt is gen[1]^j[1] ... gen[r]^j[r] */
     124             :   GEN cyc; /* t_VECSMALL of elementary divisors */
     125             : } GROUP_t;
     126             : 
     127             : static int
     128       70203 : NextElt(GROUP_t *G)
     129             : {
     130       70203 :   long i = 1;
     131       70203 :   if (G->r == 0) return 0; /* no more elt */
     132      141008 :   while (++G->j[i] == G->cyc[i]) /* from 0 to cyc[i]-1 */
     133             :   {
     134        1533 :     G->j[i] = 0;
     135        1533 :     if (++i > G->r) return 0; /* no more elt */
     136             :   }
     137       69405 :   return i; /* we have multiplied by gen[i] */
     138             : }
     139             : 
     140             : /* Compute all the elements of a group given by its SNF */
     141             : static GEN
     142        1043 : EltsOfGroup(long order, GEN cyc)
     143             : {
     144             :   long i;
     145             :   GEN rep;
     146             :   GROUP_t G;
     147             : 
     148        1043 :   G.cyc = gtovecsmall(cyc);
     149        1043 :   G.r = lg(cyc)-1;
     150        1043 :   G.j = zero_zv(G.r);
     151             : 
     152        1043 :   rep = cgetg(order + 1, t_VEC);
     153        1043 :   gel(rep,order) = vecsmall_to_col(G.j);
     154             : 
     155        7371 :   for  (i = 1; i < order; i++)
     156             :   {
     157        6328 :     (void)NextElt(&G);
     158        6328 :     gel(rep,i) = vecsmall_to_col(G.j);
     159             :   }
     160        1043 :   return rep;
     161             : }
     162             : 
     163             : /* enumerate all group elements */
     164             : GEN
     165         371 : cyc2elts(GEN cyc)
     166             : {
     167             :   long i, n;
     168             :   GEN rep;
     169             :   GROUP_t G;
     170             : 
     171         371 :   G.cyc = typ(cyc)==t_VECSMALL? cyc: gtovecsmall(cyc);
     172         371 :   n = zv_prod(G.cyc);
     173         371 :   G.r = lg(cyc)-1;
     174         371 :   G.j = zero_zv(G.r);
     175             : 
     176         371 :   rep = cgetg(n+1, t_VEC);
     177         371 :   gel(rep,n) = leafcopy(G.j); /* trivial elt comes last */
     178        1239 :   for  (i = 1; i < n; i++)
     179             :   {
     180         868 :     (void)NextElt(&G);
     181         868 :     gel(rep,i) = leafcopy(G.j);
     182             :   }
     183         371 :   return rep;
     184             : }
     185             : 
     186             : /* Let Qt as given by InitQuotient, compute a system of
     187             :    representatives of the quotient */
     188             : static GEN
     189         665 : ComputeLift(GEN Qt)
     190             : {
     191         665 :   GEN e, U = gel(Qt,3);
     192         665 :   long i, h = itos(gel(Qt,1));
     193             : 
     194         665 :   e = EltsOfGroup(h, gel(Qt,2));
     195         665 :   if (!RgM_isidentity(U))
     196             :   {
     197          28 :     GEN Ui = ZM_inv(U, gen_1);
     198          28 :     for (i = 1; i <= h; i++) gel(e,i) = ZM_ZC_mul(Ui, gel(e,i));
     199             :   }
     200         665 :   return e;
     201             : }
     202             : 
     203             : /* nchi: a character given by a vector [d, (c_i)], e.g. from char_normalize
     204             :  * such that chi(x) = e((c . log(x)) / d) where log(x) on bnr.gen */
     205             : static GEN
     206        2828 : get_Char(GEN nchi, long prec)
     207        2828 : { return mkvec2(nchi, rootsof1_cx(gel(nchi,1), prec)); }
     208             : 
     209             : /* prime divisors of conductor */
     210             : static GEN
     211         336 : divcond(GEN bnr) {GEN bid = bnr_get_bid(bnr); return gel(bid_get_fact(bid),1);}
     212             : 
     213             : /* vector of prime ideals dividing bnr but not bnrc */
     214             : static GEN
     215         210 : get_prdiff(GEN bnr, GEN condc)
     216             : {
     217         210 :   GEN prdiff, M = gel(condc,1), D = divcond(bnr), nf = bnr_get_nf(bnr);
     218         210 :   long nd, i, l  = lg(D);
     219         210 :   prdiff = cgetg(l, t_COL);
     220         616 :   for (nd=1, i=1; i < l; i++)
     221         406 :     if (!idealval(nf, M, gel(D,i))) gel(prdiff,nd++) = gel(D,i);
     222         210 :   setlg(prdiff, nd); return prdiff;
     223             : }
     224             : 
     225             : #define ch_C(x)    gel(x,1)
     226             : #define ch_bnr(x)  gel(x,2)
     227             : #define ch_4(x)    gel(x,3)
     228             : #define ch_CHI(x)  gel(x,4)
     229             : #define ch_diff(x) gel(x,5)
     230             : #define ch_cond(x) gel(x,6)
     231             : #define ch_CHI0(x) gel(x,7)
     232             : #define ch_comp(x) gel(x,8)
     233             : static long
     234        3647 : ch_deg(GEN dtcr) { return chi_get_deg(ch_CHI(dtcr)); }
     235             : 
     236             : static GEN
     237        1078 : GetDeg(GEN dataCR)
     238             : {
     239        1078 :   long i, l = lg(dataCR);
     240        1078 :   GEN degs = cgetg(l, t_VECSMALL);
     241        1078 :   for (i = 1; i < l; i++) degs[i] = eulerphiu(ch_deg(gel(dataCR,i)));
     242        1078 :   return degs;
     243             : }
     244             : 
     245             : /********************************************************************/
     246             : /*                    1rst part: find the field K                   */
     247             : /********************************************************************/
     248             : static GEN AllStark(GEN data, GEN nf, long flag, long prec);
     249             : 
     250             : /* Columns of C [HNF] give the generators of a subgroup of the finite abelian
     251             :  * group A [ in terms of implicit generators ], compute data to work in A/C:
     252             :  * 1) order
     253             :  * 2) structure
     254             :  * 3) the matrix A ->> A/C
     255             :  * 4) the subgroup C */
     256             : static GEN
     257        2058 : InitQuotient(GEN C)
     258             : {
     259        2058 :   GEN U, D = ZM_snfall_i(C, &U, NULL, 1), h = ZV_prod(D);
     260        2058 :   return mkvec5(h, D, U, C, cyc_normalize(D));
     261             : }
     262             : 
     263             : /* lift chi character on A/C [Qt from InitQuotient] to character on A [cyc]*/
     264             : static GEN
     265        2317 : LiftChar(GEN Qt, GEN cyc, GEN chi)
     266             : {
     267        2317 :   GEN ncyc = gel(Qt,5), U = gel(Qt,3);
     268        2317 :   GEN nchi = char_normalize(chi, ncyc);
     269        2317 :   GEN c = ZV_ZM_mul(gel(nchi,2), U), d = gel(nchi,1);
     270        2317 :   return char_denormalize(cyc, d, c);
     271             : }
     272             : 
     273             : /* Let s: A -> B given by P, and let cycA, cycB be the cyclic structure of
     274             :  * A and B, compute the kernel of s. */
     275             : static GEN
     276         420 : ComputeKernel0(GEN P, GEN cycA, GEN cycB)
     277             : {
     278         420 :   pari_sp av = avma;
     279         420 :   long nbA = lg(cycA)-1, rk;
     280         420 :   GEN U, DB = diagonal_shallow(cycB);
     281             : 
     282         420 :   rk = nbA + lg(cycB) - lg(ZM_hnfall_i(shallowconcat(P, DB), &U, 1));
     283         420 :   U = matslice(U, 1,nbA, 1,rk);
     284         420 :   return gerepileupto(av, ZM_hnfmodid(U, cycA));
     285             : }
     286             : 
     287             : /* Let m and n be two moduli such that n|m and let C be a congruence
     288             :    group modulo n, compute the corresponding congruence group modulo m
     289             :    ie the kernel of the map Clk(m) ->> Clk(n)/C */
     290             : static GEN
     291         420 : ComputeKernel(GEN bnrm, GEN bnrn, GEN dtQ)
     292             : {
     293         420 :   pari_sp av = avma;
     294         420 :   GEN P = ZM_mul(gel(dtQ,3), bnrsurjection(bnrm, bnrn));
     295         420 :   return gerepileupto(av, ComputeKernel0(P, bnr_get_cyc(bnrm), gel(dtQ,2)));
     296             : }
     297             : 
     298             : static long
     299        1295 : cyc_is_cyclic(GEN cyc) { return lg(cyc) <= 2 || equali1(gel(cyc,2)); }
     300             : 
     301             : /* Let H be a subgroup of cl(bnr)/sugbroup, return 1 if
     302             :    cl(bnr)/subgoup/H is cyclic and the signature of the
     303             :    corresponding field is equal to sig and no finite prime
     304             :    dividing cond(bnr) is totally split in this field. Return 0
     305             :    otherwise. */
     306             : static long
     307         644 : IsGoodSubgroup(GEN H, GEN bnr, GEN map)
     308             : {
     309         644 :   pari_sp av = avma;
     310             :   GEN mod, modH, p1, p2, U, P, PH, bnrH, iH, qH;
     311             :   long j;
     312             : 
     313         644 :   p1 = InitQuotient(H);
     314             :   /* quotient is non cyclic */
     315         644 :   if (!cyc_is_cyclic(gel(p1,2))) { avma = av; return 0; }
     316             : 
     317         294 :   p2 = ZM_hnfall_i(shallowconcat(map,H), &U, 0);
     318         294 :   setlg(U, lg(H));
     319         294 :   for (j = 1; j < lg(U); j++) setlg(gel(U,j), lg(H));
     320         294 :   p1 = ZM_hnfmodid(U, bnr_get_cyc(bnr)); /* H as a subgroup of bnr */
     321         294 :   modH = bnrconductor_i(bnr, p1, 0);
     322         294 :   mod  = bnr_get_mod(bnr);
     323             : 
     324             :   /* is the signature correct? */
     325         294 :   if (!gequal(gel(modH,2), gel(mod,2))) { avma = av; return 0; }
     326             : 
     327             :   /* finite part are the same: OK */
     328         203 :   if (gequal(gel(modH,1), gel(mod,1))) { avma = av; return 1; }
     329             : 
     330             :   /* need to check the splitting of primes dividing mod but not modH */
     331          63 :   bnrH = Buchray(bnr, modH, nf_INIT);
     332          63 :   P = divcond(bnr);
     333          63 :   PH = divcond(bnrH);
     334          63 :   p2 = ZM_mul(bnrsurjection(bnr, bnrH), p1);
     335             :   /* H as a subgroup of bnrH */
     336          63 :   iH = ZM_hnfmodid(p2,  bnr_get_cyc(bnrH));
     337          63 :   qH = InitQuotient(iH);
     338         203 :   for (j = 1; j < lg(P); j++)
     339             :   {
     340         161 :     GEN pr = gel(P, j), e;
     341             :     /* if pr divides modH, it is ramified, so it's good */
     342         161 :     if (tablesearch(PH, pr, cmp_prime_ideal)) continue;
     343             :     /* inertia degree of pr in bnr(modH)/H is charorder(e, cycH) */
     344          56 :     e = ZM_ZC_mul(gel(qH,3), isprincipalray(bnrH, pr));
     345          56 :     e = vecmodii(e, gel(qH,2));
     346          56 :     if (ZV_equal0(e)) { avma = av; return 0; } /* f = 1 */
     347             :   }
     348          42 :   avma = av; return 1;
     349             : }
     350             : 
     351             : /* compute the list of characters to consider for AllStark and
     352             :    initialize precision-independent data to compute with them */
     353             : static GEN
     354         322 : get_listCR(GEN bnr, GEN dtQ)
     355             : {
     356             :   GEN listCR, vecchi, Mr;
     357             :   long hD, h, nc, i, tnc;
     358             :   hashtable *S;
     359             : 
     360         322 :   Mr = bnr_get_cyc(bnr);
     361         322 :   hD = itos(gel(dtQ,1));
     362         322 :   h  = hD >> 1;
     363             : 
     364         322 :   listCR = cgetg(h+1, t_VEC); /* non-conjugate chars */
     365         322 :   nc = tnc = 1;
     366         322 :   vecchi = EltsOfGroup(hD, gel(dtQ,2));
     367         322 :   S = hash_create(h, (ulong(*)(void*))&hash_GEN,
     368             :                      (int(*)(void*,void*))&ZV_equal, 1);
     369        2289 :   for (i = 1; tnc <= h; i++)
     370             :   { /* lift a character of D in Clk(m) */
     371        1967 :     GEN cond, lchi = LiftChar(dtQ, Mr, gel(vecchi,i));
     372        1967 :     if (hash_search(S, lchi)) continue;
     373        1932 :     cond = bnrconductorofchar(bnr, lchi);
     374        1932 :     if (gequal0(gel(cond,2))) continue;
     375             :     /* the infinite part of chi is non trivial */
     376        1050 :     gel(listCR,nc++) = mkvec2(lchi, cond);
     377             : 
     378             :     /* if chi is not real, add its conjugate character to S */
     379        1050 :     if (absequaliu(charorder(Mr,lchi), 2)) tnc++;
     380             :     else
     381             :     {
     382         658 :       hash_insert(S, charconj(Mr, lchi), (void*)1);
     383         658 :       tnc+=2;
     384             :     }
     385             :   }
     386         322 :   setlg(listCR, nc); return listCR;
     387             : }
     388             : 
     389             : static GEN InitChar(GEN bnr, GEN listCR, long prec);
     390             : 
     391             : /* Given a conductor and a subgroups, return the corresponding
     392             :    complexity and precision required using quickpol. Fill data[5] with
     393             :    listCR */
     394             : static long
     395         322 : CplxModulus(GEN data, long *newprec)
     396             : {
     397         322 :   long pr, ex, dprec = DEFAULTPREC;
     398             :   pari_sp av;
     399         322 :   GEN pol, listCR, cpl, bnr = gel(data,1), nf = checknf(bnr);
     400             : 
     401         322 :   listCR = get_listCR(bnr, gel(data,3));
     402         322 :   for (av = avma;; avma = av)
     403             :   {
     404         322 :     gel(data,5) = InitChar(bnr, listCR, dprec);
     405         322 :     pol = AllStark(data, nf, -1, dprec);
     406         322 :     pr = nbits2extraprec( gexpo(pol) );
     407         322 :     if (pr < 0) pr = 0;
     408         322 :     dprec = maxss(dprec, pr) + EXTRA_PREC;
     409         322 :     if (!gequal0(leading_coeff(pol)))
     410             :     {
     411         322 :       cpl = RgX_fpnorml2(pol, DEFAULTPREC);
     412         322 :       if (!gequal0(cpl)) break;
     413             :     }
     414           0 :     if (DEBUGLEVEL>1) pari_warn(warnprec, "CplxModulus", dprec);
     415           0 :   }
     416         322 :   ex = gexpo(cpl); avma = av;
     417         322 :   if (DEBUGLEVEL>1) err_printf("cpl = 2^%ld\n", ex);
     418             : 
     419         322 :   gel(data,5) = listCR;
     420         322 :   *newprec = dprec; return ex;
     421             : }
     422             : 
     423             : /* return A \cap B in abelian group defined by cyc. NULL = whole group */
     424             : static GEN
     425         693 : subgp_intersect(GEN cyc, GEN A, GEN B)
     426             : {
     427             :   GEN H, U;
     428             :   long k, lH;
     429         693 :   if (!A) return B;
     430         301 :   if (!B) return A;
     431         301 :   H = ZM_hnfall_i(shallowconcat(A,B), &U, 1);
     432         301 :   setlg(U, lg(A)); lH = lg(H);
     433         301 :   for (k = 1; k < lg(U); k++) setlg(gel(U,k), lH);
     434         301 :   return ZM_hnfmodid(ZM_mul(A,U), cyc);
     435             : }
     436             : 
     437             :  /* Let f be a conductor without infinite part and let C be a
     438             :    congruence group modulo f, compute (m,D) such that D is a
     439             :    congruence group of conductor m where m is a multiple of f
     440             :    divisible by all the infinite places but one, D is a subgroup of
     441             :    index 2 of Im(C) in Clk(m), and m is such that the intersection
     442             :    of the subgroups H of Clk(m)/D such that the quotient is
     443             :    cyclic and no prime divding m, but the one infinite prime, is
     444             :    totally split in the extension corresponding to H is trivial.
     445             :    Return bnr(m), D, the quotient Ck(m)/D and Clk(m)/Im(C) */
     446             : static GEN
     447         322 : FindModulus(GEN bnr, GEN dtQ, long *newprec)
     448             : {
     449         322 :   const long limnorm = 400;
     450             :   long n, i, narch, maxnorm, minnorm, N;
     451         322 :   long first = 1, pr, rb, oldcpl = -1, iscyc;
     452         322 :   pari_sp av = avma;
     453         322 :   GEN bnf, nf, f, arch, m, rep = NULL;
     454             : 
     455         322 :   bnf = bnr_get_bnf(bnr);
     456         322 :   nf  = bnf_get_nf(bnf);
     457         322 :   N   = nf_get_degree(nf);
     458         322 :   f   = gel(bnr_get_mod(bnr), 1);
     459             : 
     460             :   /* if cpl < rb, it is not necessary to try another modulus */
     461         322 :   rb = expi( powii(mulii(nf_get_disc(nf), ZM_det_triangular(f)), gmul2n(bnr_get_no(bnr), 3)) );
     462             : 
     463             :   /* Initialization of the possible infinite part */
     464         322 :   arch = const_vec(N, gen_1);
     465             : 
     466             :   /* narch = (N == 2)? 1: N; -- if N=2, only one case is necessary */
     467         322 :   narch = N;
     468         322 :   m = mkvec2(NULL, arch);
     469             : 
     470             :   /* go from minnorm up to maxnorm. If necessary, increase these values.
     471             :    * If we cannot find a suitable conductor of norm < limnorm, stop */
     472         322 :   maxnorm = 50;
     473         322 :   minnorm = 1;
     474             : 
     475             :   /* if the extension is cyclic then we _must_ find a suitable conductor */
     476         322 :   iscyc = cyc_is_cyclic(gel(dtQ,2));
     477             : 
     478         322 :   if (DEBUGLEVEL>1)
     479           0 :     err_printf("Looking for a modulus of norm: ");
     480             : 
     481             :   for(;;)
     482             :   {
     483         322 :     GEN listid = ideallist0(nf, maxnorm, 4+8); /* ideals of norm <= maxnorm */
     484         322 :     pari_sp av1 = avma;
     485        1127 :     for (n = minnorm; n <= maxnorm; n++, avma = av1)
     486             :     {
     487        1127 :       GEN idnormn = gel(listid,n);
     488        1127 :       long nbidnn  = lg(idnormn) - 1;
     489        1127 :       if (DEBUGLEVEL>1) err_printf(" %ld", n);
     490        1813 :       for (i = 1; i <= nbidnn; i++)
     491             :       { /* finite part of the conductor */
     492             :         long s;
     493             : 
     494        1008 :         gel(m,1) = idealmul(nf, f, gel(idnormn,i));
     495        2688 :         for (s = 1; s <= narch; s++)
     496             :         { /* infinite part */
     497             :           GEN candD, ImC, bnrm;
     498             :           long nbcand, c;
     499        2002 :           gel(arch,N+1-s) = gen_0;
     500             : 
     501             :           /* compute Clk(m), check if m is a conductor */
     502        2002 :           bnrm = Buchray(bnf, m, nf_INIT);
     503        2002 :           c = bnrisconductor(bnrm, NULL);
     504        2002 :           gel(arch,N+1-s) = gen_1;
     505        2002 :           if (!c) continue;
     506             : 
     507             :           /* compute Im(C) in Clk(m)... */
     508         420 :           ImC = ComputeKernel(bnrm, bnr, dtQ);
     509             : 
     510             :           /* ... and its subgroups of index 2 with conductor m */
     511         420 :           candD = subgrouplist_cond_sub(bnrm, ImC, mkvec(gen_2));
     512         420 :           nbcand = lg(candD) - 1;
     513         427 :           for (c = 1; c <= nbcand; c++)
     514             :           {
     515         329 :             GEN D  = gel(candD,c); /* check if the conductor is suitable */
     516             :             long cpl;
     517         329 :             GEN p1 = InitQuotient(D), p2;
     518         329 :             GEN ord = gel(p1,1), cyc = gel(p1,2), map = gel(p1,3);
     519             : 
     520         329 :             if (!cyc_is_cyclic(cyc)) /* cyclic => suitable, else test */
     521             :             {
     522          84 :               GEN lH = subgrouplist(cyc, NULL), IK = NULL;
     523          84 :               long j, ok = 0;
     524         700 :               for (j = 1; j < lg(lH); j++)
     525             :               {
     526         693 :                 GEN H = gel(lH, j), IH = subgp_intersect(cyc, IK, H);
     527             :                 /* if H > IK, no need to test H */
     528         693 :                 if (IK && gidentical(IH, IK)) continue;
     529         644 :                 if (IsGoodSubgroup(H, bnrm, map))
     530             :                 {
     531         182 :                   IK = IH;
     532         182 :                   if (equalii(ord, ZM_det_triangular(IK))) { ok = 1; break; }
     533             :                 }
     534             :               }
     535          84 :               if (!ok) continue;
     536             :             }
     537             : 
     538         322 :             p2 = cgetg(6, t_VEC); /* p2[5] filled in CplxModulus */
     539         322 :             gel(p2,1) = bnrm;
     540         322 :             gel(p2,2) = D;
     541         322 :             gel(p2,3) = InitQuotient(D);
     542         322 :             gel(p2,4) = InitQuotient(ImC);
     543         322 :             if (DEBUGLEVEL>1)
     544           0 :               err_printf("\nTrying modulus = %Ps and subgroup = %Ps\n",
     545             :                          bnr_get_mod(bnrm), D);
     546         322 :             cpl = CplxModulus(p2, &pr);
     547         322 :             if (oldcpl < 0 || cpl < oldcpl)
     548             :             {
     549         322 :               *newprec = pr;
     550         322 :               if (rep) gunclone(rep);
     551         322 :               rep    = gclone(p2);
     552         322 :               oldcpl = cpl;
     553             :             }
     554         322 :             if (oldcpl < rb) goto END; /* OK */
     555             : 
     556           0 :             if (DEBUGLEVEL>1) err_printf("Trying to find another modulus...");
     557           0 :             first = 0;
     558             :           }
     559             :         }
     560         686 :         if (!first) goto END; /* OK */
     561             :       }
     562             :     }
     563             :     /* if necessary compute more ideals */
     564           0 :     minnorm = maxnorm;
     565           0 :     maxnorm <<= 1;
     566           0 :     if (!iscyc && maxnorm > limnorm) return NULL;
     567             : 
     568           0 :   }
     569             : END:
     570         322 :   if (DEBUGLEVEL>1)
     571           0 :     err_printf("No, we're done!\nModulus = %Ps and subgroup = %Ps\n",
     572           0 :                bnr_get_mod(gel(rep,1)), gel(rep,2));
     573         322 :   gel(rep,5) = InitChar(gel(rep,1), gel(rep,5), *newprec);
     574         322 :   return gerepilecopy(av, rep);
     575             : }
     576             : 
     577             : /********************************************************************/
     578             : /*                      2nd part: compute W(X)                      */
     579             : /********************************************************************/
     580             : 
     581             : /* find ilambda s.t. Diff*f*ilambda integral and coprime to f
     582             :    and ilambda >> 0 at foo, fa = factorization of f */
     583             : static GEN
     584         798 : get_ilambda(GEN nf, GEN fa, GEN foo)
     585             : {
     586         798 :   GEN x, w, E2, P = gel(fa,1), E = gel(fa,2), D = nf_get_diff(nf);
     587         798 :   long i, l = lg(P);
     588         798 :   if (l == 1) return gen_1;
     589         665 :   w = cgetg(l, t_VEC);
     590         665 :   E2 = cgetg(l, t_COL);
     591        1463 :   for (i = 1; i < l; i++)
     592             :   {
     593         798 :     GEN pr = gel(P,i), t = pr_get_tau(pr);
     594         798 :     long e = itou(gel(E,i)), v = idealval(nf, D, pr);
     595         798 :     if (v) { D = idealdivpowprime(nf, D, pr, utoipos(v)); e += v; }
     596         798 :     gel(E2,i) = stoi(e+1);
     597         798 :     if (typ(t) == t_MAT) t = gel(t,1);
     598         798 :     gel(w,i) = gdiv(nfpow(nf, t, stoi(e)), powiu(pr_get_p(pr),e));
     599             :   }
     600         665 :   x = mkmat2(P, E2);
     601         665 :   return idealchinese(nf, mkvec2(x, foo), w);
     602             : }
     603             : /* compute the list of W(chi) such that Ld(s,chi) = W(chi) Ld(1 - s, chi*),
     604             :  * for all chi in LCHI. All chi have the same conductor (= cond(bnr)).
     605             :  * if check == 0 do not check the result */
     606             : static GEN
     607         994 : ArtinNumber(GEN bnr, GEN LCHI, long check, long prec)
     608             : {
     609         994 :   long ic, i, j, nz, nChar = lg(LCHI)-1;
     610         994 :   pari_sp av = avma, av2;
     611             :   GEN sqrtnc, cond, condZ, cond0, cond1, nf, T;
     612             :   GEN cyc, vN, vB, diff, vt, idh, zid, gen, z, nchi;
     613             :   GEN indW, W, classe, s0, s, den, ilambda, sarch;
     614             :   CHI_t **lC;
     615             :   GROUP_t G;
     616             : 
     617         994 :   lC = (CHI_t**)new_chunk(nChar + 1);
     618         994 :   indW = cgetg(nChar + 1, t_VECSMALL);
     619         994 :   W = cgetg(nChar + 1, t_VEC);
     620        3472 :   for (ic = 0, i = 1; i <= nChar; i++)
     621             :   {
     622        2478 :     GEN CHI = gel(LCHI,i);
     623        2478 :     if (chi_get_deg(CHI) <= 2) { gel(W,i) = gen_1; continue; }
     624        1533 :     ic++; indW[ic] = i;
     625        1533 :     lC[ic] = (CHI_t*)new_chunk(sizeof(CHI_t));
     626        1533 :     init_CHI_C(lC[ic], CHI);
     627             :   }
     628         994 :   if (!ic) return W;
     629         798 :   nChar = ic;
     630             : 
     631         798 :   nf    = bnr_get_nf(bnr);
     632         798 :   diff  = nf_get_diff(nf);
     633         798 :   T     = nf_get_Tr(nf);
     634         798 :   cond  = bnr_get_mod(bnr);
     635         798 :   cond0 = gel(cond,1); condZ = gcoeff(cond0,1,1);
     636         798 :   cond1 = gel(cond,2);
     637             : 
     638         798 :   sqrtnc = gsqrt(idealnorm(nf, cond0), prec);
     639         798 :   ilambda = get_ilambda(nf, bid_get_fact(bnr_get_bid(bnr)), cond1);
     640         798 :   idh = idealmul(nf, ilambda, idealmul(nf, diff, cond0)); /* integral */
     641         798 :   ilambda = Q_remove_denom(ilambda, &den);
     642         798 :   z = den? rootsof1_cx(den, prec): NULL;
     643             : 
     644             :   /* compute a system of generators of (Ok/cond)^*, we'll make them
     645             :    * cond1-positive in the main loop */
     646         798 :   zid = Idealstar(nf, cond0, nf_GEN);
     647         798 :   cyc = abgrp_get_cyc(zid);
     648         798 :   gen = abgrp_get_gen(zid);
     649         798 :   nz = lg(gen) - 1;
     650         798 :   sarch = nfarchstar(nf, cond0, vec01_to_indices(cond1));
     651             : 
     652         798 :   nchi = cgetg(nChar+1, t_VEC);
     653         798 :   for (ic = 1; ic <= nChar; ic++) gel(nchi,ic) = cgetg(nz + 1, t_VECSMALL);
     654        1687 :   for (i = 1; i <= nz; i++)
     655             :   {
     656         889 :     if (is_bigint(gel(cyc,i)))
     657           0 :       pari_err_OVERFLOW("ArtinNumber [conductor too large]");
     658         889 :     gel(gen,i) = set_sign_mod_divisor(nf, NULL, gel(gen,i), sarch);
     659         889 :     classe = isprincipalray(bnr, gel(gen,i));
     660        2681 :     for (ic = 1; ic <= nChar; ic++) {
     661        1792 :       GEN n = gel(nchi,ic);
     662        1792 :       n[i] = CHI_eval_n(lC[ic], classe);
     663             :     }
     664             :   }
     665             : 
     666             :   /* Sum chi(beta) * exp(2i * Pi * Tr(beta * ilambda) where beta
     667             :      runs through the classes of (Ok/cond0)^* and beta cond1-positive */
     668         798 :   vt = gel(T,1); /* ( Tr(w_i) )_i */
     669         798 :   if (typ(ilambda) == t_COL)
     670         665 :     vt = ZV_ZM_mul(vt, zk_multable(nf, ilambda));
     671             :   else
     672         133 :     vt = ZC_Z_mul(vt, ilambda);
     673             :   /*vt = den . (Tr(w_i * ilambda))_i */
     674         798 :   G.cyc = gtovecsmall(cyc);
     675         798 :   G.r = nz;
     676         798 :   G.j = zero_zv(nz);
     677         798 :   vN = zero_Flm_copy(nz, nChar);
     678             : 
     679         798 :   av2 = avma;
     680         798 :   vB = const_vec(nz, gen_1);
     681         798 :   s0 = z? powgi(z, modii(gel(vt,1), den)): gen_1; /* for beta = 1 */
     682         798 :   s = const_vec(nChar, s0);
     683             : 
     684       63805 :   while ( (i = NextElt(&G)) )
     685             :   {
     686       62209 :     GEN b = gel(vB,i);
     687       62209 :     b = nfmuli(nf, b, gel(gen,i));
     688       62209 :     b = typ(b) == t_COL? FpC_red(b, condZ): modii(b, condZ);
     689       62209 :     for (j=1; j<=i; j++) gel(vB,j) = b;
     690             : 
     691      256746 :     for (ic = 1; ic <= nChar; ic++)
     692             :     {
     693      194537 :       GEN v = gel(vN,ic), n = gel(nchi,ic);
     694      194537 :       v[i] = Fl_add(v[i], n[i], lC[ic]->ord);
     695      194537 :       for (j=1; j<i; j++) v[j] = v[i];
     696             :     }
     697             : 
     698       62209 :     gel(vB,i) = b = set_sign_mod_divisor(nf, NULL, b, sarch);
     699       62209 :     if (!z)
     700           0 :       s0 = gen_1;
     701             :     else
     702             :     {
     703       62209 :       b = typ(b) == t_COL? ZV_dotproduct(vt, b): mulii(gel(vt,1),b);
     704       62209 :       s0 = powgi(z, modii(b,den));
     705             :     }
     706      256746 :     for (ic = 1; ic <= nChar; ic++)
     707             :     {
     708      194537 :       GEN v = gel(vN,ic), val = lC[ic]->val[ v[i] ];
     709      194537 :       gel(s,ic) = gadd(gel(s,ic), gmul(val, s0));
     710             :     }
     711             : 
     712       62209 :     if (gc_needed(av2, 1))
     713             :     {
     714           1 :       if (DEBUGMEM > 1) pari_warn(warnmem,"ArtinNumber");
     715           1 :       gerepileall(av2, 2, &s, &vB);
     716             :     }
     717             :   }
     718             : 
     719         798 :   classe = isprincipalray(bnr, idh);
     720         798 :   z = powIs(- (lg(gel(sarch,1))-1));
     721             : 
     722        2331 :   for (ic = 1; ic <= nChar; ic++)
     723             :   {
     724        1533 :     s0 = gmul(gel(s,ic), CHI_eval(lC[ic], classe));
     725        1533 :     s0 = gdiv(s0, sqrtnc);
     726        1533 :     if (check && - expo(subrs(gnorm(s0), 1)) < prec2nbits(prec) >> 1)
     727           0 :       pari_err_BUG("ArtinNumber");
     728        1533 :     gel(W, indW[ic]) = gmul(s0, z);
     729             :   }
     730         798 :   return gerepilecopy(av, W);
     731             : }
     732             : 
     733             : static GEN
     734         721 : ComputeAllArtinNumbers(GEN dataCR, GEN vChar, int check, long prec)
     735             : {
     736         721 :   long j, k, cl = lg(dataCR) - 1, J = lg(vChar)-1;
     737         721 :   GEN W = cgetg(cl+1,t_VEC), WbyCond, LCHI;
     738             : 
     739        1638 :   for (j = 1; j <= J; j++)
     740             :   {
     741         917 :     GEN LChar = gel(vChar,j), ldata = vecpermute(dataCR, LChar);
     742         917 :     GEN dtcr = gel(ldata,1), bnr = ch_bnr(dtcr);
     743         917 :     long l = lg(LChar);
     744             : 
     745         917 :     if (DEBUGLEVEL>1)
     746           0 :       err_printf("* Root Number: cond. no %ld/%ld (%ld chars)\n", j, J, l-1);
     747         917 :     LCHI = cgetg(l, t_VEC);
     748         917 :     for (k = 1; k < l; k++) gel(LCHI,k) = ch_CHI0(gel(ldata,k));
     749         917 :     WbyCond = ArtinNumber(bnr, LCHI, check, prec);
     750         917 :     for (k = 1; k < l; k++) gel(W,LChar[k]) = gel(WbyCond,k);
     751             :   }
     752         721 :   return W;
     753             : }
     754             : static GEN
     755          77 : SingleArtinNumber(GEN bnr, GEN chi, long prec)
     756          77 : { return gel(ArtinNumber(bnr, mkvec(chi), 1, prec), 1); }
     757             : 
     758             : /* compute the constant W of the functional equation of
     759             :    Lambda(chi). If flag = 1 then chi is assumed to be primitive */
     760             : GEN
     761          77 : bnrrootnumber(GEN bnr, GEN chi, long flag, long prec)
     762             : {
     763          77 :   pari_sp av = avma;
     764             :   GEN cyc;
     765             : 
     766          77 :   if (flag < 0 || flag > 1) pari_err_FLAG("bnrrootnumber");
     767          77 :   checkbnr(bnr);
     768          77 :   if (flag)
     769             :   {
     770           0 :     cyc = bnr_get_cyc(bnr);
     771           0 :     if (!char_check(cyc,chi)) pari_err_TYPE("bnrrootnumber [character]", chi);
     772             :   }
     773             :   else
     774             :   {
     775          77 :     GEN z = bnrconductor_i(bnr, chi, 2);
     776          77 :     bnr = gel(z,2);
     777          77 :     chi = gel(z,3);
     778          77 :     cyc = bnr_get_cyc(bnr);
     779             :   }
     780          77 :   chi = char_normalize(chi, cyc_normalize(cyc));
     781          77 :   chi = get_Char(chi, prec);
     782          77 :   return gerepilecopy(av, SingleArtinNumber(bnr, chi, prec));
     783             : }
     784             : 
     785             : /********************************************************************/
     786             : /*               3rd part: initialize the characters                */
     787             : /********************************************************************/
     788             : 
     789             : /* Let chi be a character, A(chi) corresponding to the primes dividing diff
     790             :    at s = flag. If s = 0, returns [r, A] where r is the order of vanishing
     791             :    at s = 0 corresponding to diff. No GC */
     792             : static GEN
     793        2163 : ComputeAChi(GEN dtcr, long *r, long flag, long prec)
     794             : {
     795        2163 :   GEN A, diff = ch_diff(dtcr), bnrc = ch_bnr(dtcr), chi  = ch_CHI0(dtcr);
     796        2163 :   long i, l = lg(diff);
     797             : 
     798        2163 :   A = gen_1; *r = 0;
     799        2282 :   for (i = 1; i < l; i++)
     800             :   {
     801         119 :     GEN pr = gel(diff,i), B;
     802         119 :     GEN z = CharEval(chi, isprincipalray(bnrc, pr));
     803             : 
     804         119 :     if (flag)
     805           0 :       B = gsubsg(1, gdiv(z, pr_norm(pr)));
     806         119 :     else if (gequal1(z))
     807             :     {
     808          14 :       B = glog(pr_norm(pr), prec);
     809          14 :       (*r)++;
     810             :     }
     811             :     else
     812         105 :       B = gsubsg(1, z);
     813         119 :     A = gmul(A, B);
     814             :   }
     815        2163 :   return A;
     816             : }
     817             : /* simplified version of ComputeAchi: return 1 if L(0,chi) = 0 */
     818             : static int
     819        2177 : L_vanishes_at_0(GEN dtcr)
     820             : {
     821        2177 :   GEN diff = ch_diff(dtcr), bnrc = ch_bnr(dtcr), chi  = ch_CHI0(dtcr);
     822        2177 :   long i, l = lg(diff);
     823             : 
     824        2282 :   for (i = 1; i < l; i++)
     825             :   {
     826         133 :     GEN pr = gel(diff,i);
     827         133 :     if (! CharEval_n(chi, isprincipalray(bnrc, pr))) return 1;
     828             :   }
     829        2149 :   return 0;
     830             : }
     831             : 
     832             : static GEN
     833         889 : _data4(GEN arch, long r1, long r2)
     834             : {
     835         889 :   GEN z = cgetg(5, t_VECSMALL);
     836         889 :   long i, b, q = 0;
     837             : 
     838         889 :   for (i=1; i<=r1; i++) if (signe(gel(arch,i))) q++;
     839         889 :   z[1] = q; b = r1 - q;
     840         889 :   z[2] = b;
     841         889 :   z[3] = r2;
     842         889 :   z[4] = maxss(b+r2+1, r2+q);
     843         889 :   return z;
     844             : }
     845             : 
     846             : /* Given a list [chi, F = cond(chi)] of characters over Cl(bnr), compute a
     847             :    vector dataCR containing for each character:
     848             :    2: the constant C(F) [t_REAL]
     849             :    3: bnr(F)
     850             :    4: [q, r1 - q, r2, rc] where
     851             :         q = number of real places in F
     852             :         rc = max{r1 + r2 - q + 1, r2 + q}
     853             :    6: diff(chi) primes dividing m but not F
     854             :    7: finite part of F
     855             : 
     856             :    1: chi
     857             :    5: [(c_i), z, d] in bnr(m)
     858             :    8: [(c_i), z, d] in bnr(F)
     859             :    9: if NULL then does not compute (for AllStark) */
     860             : static GEN
     861         700 : InitChar(GEN bnr, GEN listCR, long prec)
     862             : {
     863         700 :   GEN bnf = checkbnf(bnr), nf = bnf_get_nf(bnf);
     864             :   GEN modul, dk, C, dataCR, chi, cond, ncyc;
     865             :   long N, r1, r2, prec2, i, j, l;
     866         700 :   pari_sp av = avma;
     867             : 
     868         700 :   modul = bnr_get_mod(bnr);
     869         700 :   dk    = nf_get_disc(nf);
     870         700 :   N     = nf_get_degree(nf);
     871         700 :   nf_get_sign(nf, &r1,&r2);
     872         700 :   prec2 = precdbl(prec) + EXTRA_PREC;
     873         700 :   C     = gmul2n(sqrtr_abs(divir(dk, powru(mppi(prec2),N))), -r2);
     874         700 :   ncyc = cyc_normalize( bnr_get_cyc(bnr) );
     875             : 
     876         700 :   dataCR = cgetg_copy(listCR, &l);
     877        3024 :   for (i = 1; i < l; i++)
     878             :   {
     879        2324 :     GEN bnrc, olddtcr, dtcr = cgetg(9, t_VEC);
     880        2324 :     gel(dataCR,i) = dtcr;
     881             : 
     882        2324 :     chi  = gmael(listCR, i, 1);
     883        2324 :     cond = gmael(listCR, i, 2);
     884             : 
     885             :     /* do we already know the invariants of chi? */
     886        2324 :     olddtcr = NULL;
     887        3178 :     for (j = 1; j < i; j++)
     888        2289 :       if (gequal(cond, gmael(listCR,j,2))) { olddtcr = gel(dataCR,j); break; }
     889             : 
     890        2324 :     if (!olddtcr)
     891             :     {
     892         889 :       ch_C(dtcr) = gmul(C, gsqrt(ZM_det_triangular(gel(cond,1)), prec2));
     893         889 :       ch_4(dtcr) = _data4(gel(cond,2),r1,r2);
     894         889 :       ch_cond(dtcr) = cond;
     895         889 :       if (gequal(cond,modul))
     896             :       {
     897         679 :         ch_bnr(dtcr) = bnr;
     898         679 :         ch_diff(dtcr) = cgetg(1, t_VEC);
     899             :       }
     900             :       else
     901             :       {
     902         210 :         ch_bnr(dtcr) = Buchray(bnf, cond, nf_INIT);
     903         210 :         ch_diff(dtcr) = get_prdiff(bnr, cond);
     904             :       }
     905             :     }
     906             :     else
     907             :     {
     908        1435 :       ch_C(dtcr) = ch_C(olddtcr);
     909        1435 :       ch_bnr(dtcr) = ch_bnr(olddtcr);
     910        1435 :       ch_4(dtcr) = ch_4(olddtcr);
     911        1435 :       ch_diff(dtcr) = ch_diff(olddtcr);
     912        1435 :       ch_cond(dtcr) = ch_cond(olddtcr);
     913             :     }
     914             : 
     915        2324 :     chi = char_normalize(chi,ncyc);
     916        2324 :     ch_CHI(dtcr) = get_Char(chi, prec2);
     917        2324 :     ch_comp(dtcr) = gen_1; /* compute this character (by default) */
     918             : 
     919        2324 :     bnrc = ch_bnr(dtcr);
     920        2324 :     if (gequal(bnr_get_mod(bnr), bnr_get_mod(bnrc)))
     921        2051 :       ch_CHI0(dtcr) = ch_CHI(dtcr);
     922             :     else
     923             :     {
     924         273 :       chi = bnrchar_primitive(bnr, chi, bnrc);
     925         273 :       ch_CHI0(dtcr) = get_Char(chi, prec2);
     926             :     }
     927             :   }
     928             : 
     929         700 :   return gerepilecopy(av, dataCR);
     930             : }
     931             : 
     932             : /* recompute dataCR with the new precision */
     933             : static GEN
     934          21 : CharNewPrec(GEN dataCR, GEN nf, long prec)
     935             : {
     936             :   GEN dk, C;
     937             :   long N, l, j, prec2;
     938             : 
     939          21 :   dk    =  nf_get_disc(nf);
     940          21 :   N     =  nf_get_degree(nf);
     941          21 :   prec2 = precdbl(prec) + EXTRA_PREC;
     942             : 
     943          21 :   C = sqrtr(divir(absi(dk), powru(mppi(prec2), N)));
     944             : 
     945          21 :   l = lg(dataCR);
     946          98 :   for (j = 1; j < l; j++)
     947             :   {
     948          77 :     GEN dtcr = gel(dataCR,j), f0 = gel(ch_cond(dtcr),1);
     949          77 :     ch_C(dtcr) = gmul(C, gsqrt(ZM_det_triangular(f0), prec2));
     950             : 
     951          77 :     gmael(ch_bnr(dtcr), 1, 7) = nf;
     952             : 
     953          77 :     ch_CHI( dtcr) = get_Char(gel(ch_CHI(dtcr), 1), prec2);
     954          77 :     ch_CHI0(dtcr) = get_Char(gel(ch_CHI0(dtcr),1), prec2);
     955             :   }
     956             : 
     957          21 :   return dataCR;
     958             : }
     959             : 
     960             : /********************************************************************/
     961             : /*             4th part: compute the coefficients an(chi)           */
     962             : /*                                                                  */
     963             : /* matan entries are arrays of ints containing the coefficients of  */
     964             : /* an(chi) as a polmod modulo polcyclo(order(chi))                     */
     965             : /********************************************************************/
     966             : 
     967             : static void
     968     1506138 : _0toCoeff(int *rep, long deg)
     969             : {
     970             :   long i;
     971     1506138 :   for (i=0; i<deg; i++) rep[i] = 0;
     972     1506138 : }
     973             : 
     974             : /* transform a polmod into Coeff */
     975             : static void
     976      383425 : Polmod2Coeff(int *rep, GEN polmod, long deg)
     977             : {
     978             :   long i;
     979      383425 :   if (typ(polmod) == t_POLMOD)
     980             :   {
     981      271399 :     GEN pol = gel(polmod,2);
     982      271399 :     long d = degpol(pol);
     983             : 
     984      271399 :     pol += 2;
     985      271399 :     for (i=0; i<=d; i++) rep[i] = itos(gel(pol,i));
     986      271399 :     for (   ; i<deg; i++) rep[i] = 0;
     987             :   }
     988             :   else
     989             :   {
     990      112026 :     rep[0] = itos(polmod);
     991      112026 :     for (i=1; i<deg; i++) rep[i] = 0;
     992             :   }
     993      383425 : }
     994             : 
     995             : /* initialize a deg * n matrix of ints */
     996             : static int**
     997        4158 : InitMatAn(long n, long deg, long flag)
     998             : {
     999             :   long i, j;
    1000        4158 :   int *a, **A = (int**)pari_malloc((n+1)*sizeof(int*));
    1001        4158 :   A[0] = NULL;
    1002     5884825 :   for (i = 1; i <= n; i++)
    1003             :   {
    1004     5880667 :     a = (int*)pari_malloc(deg*sizeof(int));
    1005     5880667 :     A[i] = a; a[0] = (i == 1 || flag);
    1006     5880667 :     for (j = 1; j < deg; j++) a[j] = 0;
    1007             :   }
    1008        4158 :   return A;
    1009             : }
    1010             : 
    1011             : static void
    1012        6545 : FreeMat(int **A, long n)
    1013             : {
    1014             :   long i;
    1015     5898811 :   for (i = 0; i <= n; i++)
    1016     5892266 :     if (A[i]) pari_free((void*)A[i]);
    1017        6545 :   pari_free((void*)A);
    1018        6545 : }
    1019             : 
    1020             : /* initialize Coeff reduction */
    1021             : static int**
    1022        2387 : InitReduction(long d, long deg)
    1023             : {
    1024             :   long j;
    1025        2387 :   pari_sp av = avma;
    1026             :   int **A;
    1027             :   GEN polmod, pol;
    1028             : 
    1029        2387 :   A   = (int**)pari_malloc(deg*sizeof(int*));
    1030        2387 :   pol = polcyclo(d, 0);
    1031        9828 :   for (j = 0; j < deg; j++)
    1032             :   {
    1033        7441 :     A[j] = (int*)pari_malloc(deg*sizeof(int));
    1034        7441 :     polmod = gmodulo(pol_xn(deg+j, 0), pol);
    1035        7441 :     Polmod2Coeff(A[j], polmod, deg);
    1036             :   }
    1037             : 
    1038        2387 :   avma = av; return A;
    1039             : }
    1040             : 
    1041             : #if 0
    1042             : void
    1043             : pan(int **an, long n, long deg)
    1044             : {
    1045             :   long i,j;
    1046             :   for (i = 1; i <= n; i++)
    1047             :   {
    1048             :     err_printf("n = %ld: ",i);
    1049             :     for (j = 0; j < deg; j++) err_printf("%d ",an[i][j]);
    1050             :     err_printf("\n");
    1051             :   }
    1052             : }
    1053             : #endif
    1054             : 
    1055             : /* returns 0 if c is zero, 1 otherwise. */
    1056             : static int
    1057     7208197 : IsZero(int* c, long deg)
    1058             : {
    1059             :   long i;
    1060    22963986 :   for (i = 0; i < deg; i++)
    1061    17633591 :     if (c[i]) return 0;
    1062     5330395 :   return 1;
    1063             : }
    1064             : 
    1065             : /* set c0 <-- c0 * c1 */
    1066             : static void
    1067     1803378 : MulCoeff(int *c0, int* c1, int** reduc, long deg)
    1068             : {
    1069             :   long i,j;
    1070             :   int c, *T;
    1071             : 
    1072     3606756 :   if (IsZero(c0,deg)) return;
    1073             : 
    1074      802035 :   T = (int*)new_chunk(2*deg);
    1075     7431979 :   for (i = 0; i < 2*deg; i++)
    1076             :   {
    1077     6629944 :     c = 0;
    1078    50013504 :     for (j = 0; j <= i; j++)
    1079    43383560 :       if (j < deg && j > i - deg) c += c0[j] * c1[i-j];
    1080     6629944 :     T[i] = c;
    1081             :   }
    1082     4117007 :   for (i = 0; i < deg; i++)
    1083             :   {
    1084     3314972 :     c = T[i];
    1085     3314972 :     for (j = 0; j < deg; j++) c += reduc[j][i] * T[deg+j];
    1086     3314972 :     c0[i] = c;
    1087             :   }
    1088             : }
    1089             : 
    1090             : /* c0 <- c0 + c1 * c2 */
    1091             : static void
    1092     5404819 : AddMulCoeff(int *c0, int *c1, int* c2, int** reduc, long deg)
    1093             : {
    1094             :   long i, j;
    1095             :   pari_sp av;
    1096             :   int c, *t;
    1097             : 
    1098     5404819 :   if (IsZero(c2,deg)) return;
    1099     1075767 :   if (!c1) /* c1 == 1 */
    1100             :   {
    1101      290451 :     for (i = 0; i < deg; i++) c0[i] += c2[i];
    1102      290451 :     return;
    1103             :   }
    1104      785316 :   av = avma;
    1105      785316 :   t = (int*)new_chunk(2*deg); /* = c1 * c2, not reduced */
    1106     5432882 :   for (i = 0; i < 2*deg; i++)
    1107             :   {
    1108     4647566 :     c = 0;
    1109    28288575 :     for (j = 0; j <= i; j++)
    1110    23641009 :       if (j < deg && j > i - deg) c += c1[j] * c2[i-j];
    1111     4647566 :     t[i] = c;
    1112             :   }
    1113     3109099 :   for (i = 0; i < deg; i++)
    1114             :   {
    1115     2323783 :     c = t[i];
    1116     2323783 :     for (j = 0; j < deg; j++) c += reduc[j][i] * t[deg+j];
    1117     2323783 :     c0[i] += c;
    1118             :   }
    1119      785316 :   avma = av;
    1120             : }
    1121             : 
    1122             : /* evaluate the Coeff. No Garbage collector */
    1123             : static GEN
    1124     3701998 : EvalCoeff(GEN z, int* c, long deg)
    1125             : {
    1126             :   long i,j;
    1127             :   GEN e, r;
    1128             : 
    1129     3701998 :   if (!c) return gen_0;
    1130             : #if 0
    1131             :   /* standard Horner */
    1132             :   e = stoi(c[deg - 1]);
    1133             :   for (i = deg - 2; i >= 0; i--)
    1134             :     e = gadd(stoi(c[i]), gmul(z, e));
    1135             : #else
    1136             :   /* specific attention to sparse polynomials */
    1137     3701998 :   e = NULL;
    1138     5220452 :   for (i = deg-1; i >=0; i=j-1)
    1139             :   {
    1140    11658972 :     for (j=i; c[j] == 0; j--)
    1141    10140518 :       if (j==0)
    1142             :       {
    1143     3048844 :         if (!e) return NULL;
    1144      293462 :         if (i!=j) z = gpowgs(z,i-j+1);
    1145      293462 :         return gmul(e,z);
    1146             :       }
    1147     1518454 :     if (e)
    1148             :     {
    1149      571838 :       r = (i==j)? z: gpowgs(z,i-j+1);
    1150      571838 :       e = gadd(gmul(e,r), stoi(c[j]));
    1151             :     }
    1152             :     else
    1153      946616 :       e = stoi(c[j]);
    1154             :   }
    1155             : #endif
    1156      653154 :   return e;
    1157             : }
    1158             : 
    1159             : /* copy the n * (m+1) array matan */
    1160             : static void
    1161      317674 : CopyCoeff(int** a, int** a2, long n, long m)
    1162             : {
    1163             :   long i,j;
    1164             : 
    1165     4636016 :   for (i = 1; i <= n; i++)
    1166             :   {
    1167     4318342 :     int *b = a[i], *b2 = a2[i];
    1168     4318342 :     for (j = 0; j < m; j++) b2[j] = b[j];
    1169             :   }
    1170      317674 : }
    1171             : 
    1172             : /* return q*p if <= n. Beware overflow */
    1173             : static long
    1174      611978 : next_pow(long q, long p, long n)
    1175             : {
    1176      611978 :   const GEN x = muluu((ulong)q, (ulong)p);
    1177      611978 :   const ulong qp = uel(x,2);
    1178      611978 :   return (lgefint(x) > 3 || qp > (ulong)n)? 0: qp;
    1179             : }
    1180             : 
    1181             : static void
    1182      317674 : an_AddMul(int **an,int **an2, long np, long n, long deg, GEN chi, int **reduc)
    1183             : {
    1184      317674 :   GEN chi2 = chi;
    1185             :   long q, qk, k;
    1186      317674 :   int *c, *c2 = (int*)new_chunk(deg);
    1187             : 
    1188      317674 :   CopyCoeff(an, an2, n/np, deg);
    1189      317674 :   for (q=np;;)
    1190             :   {
    1191      349566 :     if (gequal1(chi2)) c = NULL; else { Polmod2Coeff(c2, chi2, deg); c = c2; }
    1192     5754385 :     for(k = 1, qk = q; qk <= n; k++, qk += q)
    1193     5404819 :       AddMulCoeff(an[qk], c, an2[k], reduc, deg);
    1194      349566 :     if (! (q = next_pow(q,np, n)) ) break;
    1195             : 
    1196       31892 :     chi2 = gmul(chi2, chi);
    1197       31892 :   }
    1198      317674 : }
    1199             : 
    1200             : /* correct the coefficients an(chi) according with diff(chi) in place */
    1201             : static void
    1202        2387 : CorrectCoeff(GEN dtcr, int** an, int** reduc, long n, long deg)
    1203             : {
    1204        2387 :   pari_sp av = avma;
    1205             :   long lg, j;
    1206             :   pari_sp av1;
    1207             :   int **an2;
    1208             :   GEN bnrc, diff;
    1209             :   CHI_t C;
    1210             : 
    1211        2387 :   diff = ch_diff(dtcr); lg = lg(diff) - 1;
    1212        4774 :   if (!lg) return;
    1213             : 
    1214         175 :   if (DEBUGLEVEL>2) err_printf("diff(CHI) = %Ps", diff);
    1215         175 :   bnrc = ch_bnr(dtcr);
    1216         175 :   init_CHI_alg(&C, ch_CHI0(dtcr));
    1217             : 
    1218         175 :   an2 = InitMatAn(n, deg, 0);
    1219         175 :   av1 = avma;
    1220         385 :   for (j = 1; j <= lg; j++)
    1221             :   {
    1222         210 :     GEN pr = gel(diff,j);
    1223         210 :     long Np = upr_norm(pr);
    1224         210 :     GEN chi  = CHI_eval(&C, isprincipalray(bnrc, pr));
    1225         210 :     an_AddMul(an,an2,Np,n,deg,chi,reduc);
    1226         210 :     avma = av1;
    1227             :   }
    1228         175 :   FreeMat(an2, n); avma = av;
    1229             : }
    1230             : 
    1231             : /* compute the coefficients an in the general case */
    1232             : static int**
    1233        1596 : ComputeCoeff(GEN dtcr, LISTray *R, long n, long deg)
    1234             : {
    1235        1596 :   pari_sp av = avma, av2;
    1236             :   long i, l;
    1237             :   int **an, **reduc, **an2;
    1238             :   GEN L;
    1239             :   CHI_t C;
    1240             : 
    1241        1596 :   init_CHI_alg(&C, ch_CHI(dtcr));
    1242        1596 :   an  = InitMatAn(n, deg, 0);
    1243        1596 :   an2 = InitMatAn(n, deg, 0);
    1244        1596 :   reduc  = InitReduction(C.ord, deg);
    1245        1596 :   av2 = avma;
    1246             : 
    1247        1596 :   L = R->L1; l = lg(L);
    1248      319060 :   for (i=1; i<l; i++, avma = av2)
    1249             :   {
    1250      317464 :     long np = L[i];
    1251      317464 :     GEN chi  = CHI_eval(&C, gel(R->L1ray,i));
    1252      317464 :     an_AddMul(an,an2,np,n,deg,chi,reduc);
    1253             :   }
    1254        1596 :   FreeMat(an2, n);
    1255             : 
    1256        1596 :   CorrectCoeff(dtcr, an, reduc, n, deg);
    1257        1596 :   FreeMat(reduc, deg-1);
    1258        1596 :   avma = av; return an;
    1259             : }
    1260             : 
    1261             : /********************************************************************/
    1262             : /*              5th part: compute L-functions at s=1                */
    1263             : /********************************************************************/
    1264             : static void
    1265         420 : deg11(LISTray *R, long p, GEN bnr, GEN pr) {
    1266         420 :   GEN z = isprincipalray(bnr, pr);
    1267         420 :   vecsmalltrunc_append(R->L1, p);
    1268         420 :   vectrunc_append(R->L1ray, z);
    1269         420 : }
    1270             : static void
    1271       32199 : deg12(LISTray *R, long p, GEN bnr, GEN Lpr) {
    1272       32199 :   GEN z = isprincipalray(bnr, gel(Lpr,1));
    1273       32199 :   vecsmalltrunc_append(R->L11, p);
    1274       32199 :   vectrunc_append(R->L11ray, z);
    1275       32199 : }
    1276             : static void
    1277          70 : deg0(LISTray *R, long p) { vecsmalltrunc_append(R->L0, p); }
    1278             : static void
    1279       34116 : deg2(LISTray *R, long p) { vecsmalltrunc_append(R->L2, p); }
    1280             : 
    1281             : static void
    1282         217 : InitPrimesQuad(GEN bnr, ulong N0, LISTray *R)
    1283             : {
    1284         217 :   pari_sp av = avma;
    1285         217 :   GEN bnf = bnr_get_bnf(bnr), cond = gel(bnr_get_mod(bnr), 1);
    1286         217 :   long p,i,l, condZ = itos(gcoeff(cond,1,1)), contZ = itos(content(cond));
    1287         217 :   GEN prime, Lpr, nf = bnf_get_nf(bnf), dk = nf_get_disc(nf);
    1288             :   forprime_t T;
    1289             : 
    1290         217 :   l = 1 + primepi_upper_bound(N0);
    1291         217 :   R->L0 = vecsmalltrunc_init(l);
    1292         217 :   R->L2 = vecsmalltrunc_init(l); R->condZ = condZ;
    1293         217 :   R->L1 = vecsmalltrunc_init(l); R->L1ray = vectrunc_init(l);
    1294         217 :   R->L11= vecsmalltrunc_init(l); R->L11ray= vectrunc_init(l);
    1295         217 :   prime = utoipos(2);
    1296         217 :   u_forprime_init(&T, 2, N0);
    1297       67239 :   while ( (p = u_forprime_next(&T)) )
    1298             :   {
    1299       66805 :     prime[2] = p;
    1300       66805 :     switch (kroiu(dk, p))
    1301             :     {
    1302             :     case -1: /* inert */
    1303       34137 :       if (condZ % p == 0) deg0(R,p); else deg2(R,p);
    1304       34137 :       break;
    1305             :     case 1: /* split */
    1306       32374 :       Lpr = idealprimedec(nf, prime);
    1307       32374 :       if      (condZ % p != 0) deg12(R, p, bnr, Lpr);
    1308         175 :       else if (contZ % p == 0) deg0(R,p);
    1309             :       else
    1310             :       {
    1311         175 :         GEN pr = idealval(nf, cond, gel(Lpr,1))? gel(Lpr,2): gel(Lpr,1);
    1312         175 :         deg11(R, p, bnr, pr);
    1313             :       }
    1314       32374 :       break;
    1315             :     default: /* ramified */
    1316         294 :       if (condZ % p == 0) deg0(R,p);
    1317             :       else {
    1318         245 :         GEN pr = gel(idealprimedec(nf,prime),1);
    1319         245 :         deg11(R, p, bnr, pr);
    1320             :       }
    1321         294 :       break;
    1322             :     }
    1323             :   }
    1324             :   /* precompute isprincipalray(x), x in Z */
    1325         217 :   R->rayZ = cgetg(condZ, t_VEC);
    1326        3388 :   for (i=1; i<condZ; i++)
    1327        3171 :     gel(R->rayZ,i) = (ugcd(i,condZ) == 1)? isprincipalray(bnr, utoipos(i)): gen_0;
    1328         217 :   gerepileall(av, 7, &(R->L0), &(R->L2), &(R->rayZ),
    1329             :               &(R->L1), &(R->L1ray), &(R->L11), &(R->L11ray) );
    1330         217 : }
    1331             : 
    1332             : static void
    1333         504 : InitPrimes(GEN bnr, ulong N0, LISTray *R)
    1334             : {
    1335         504 :   GEN bnf = bnr_get_bnf(bnr), cond = gel(bnr_get_mod(bnr), 1);
    1336         504 :   long p,j,k,l, condZ = itos(gcoeff(cond,1,1)), N = lg(cond)-1;
    1337         504 :   GEN tmpray, tabpr, prime, BOUND, nf = bnf_get_nf(bnf);
    1338             :   forprime_t T;
    1339             : 
    1340         504 :   R->condZ = condZ; l = primepi_upper_bound(N0) * N;
    1341         504 :   tmpray = cgetg(N+1, t_VEC);
    1342         504 :   R->L1 = vecsmalltrunc_init(l);
    1343         504 :   R->L1ray = vectrunc_init(l);
    1344         504 :   u_forprime_init(&T, 2, N0);
    1345         504 :   prime = utoipos(2);
    1346         504 :   BOUND = utoi(N0);
    1347      117411 :   while ( (p = u_forprime_next(&T)) )
    1348             :   {
    1349      116403 :     pari_sp av = avma;
    1350      116403 :     prime[2] = p;
    1351      116403 :     if (DEBUGLEVEL>1 && (p & 2047) == 1) err_printf("%ld ", p);
    1352      116403 :     tabpr = idealprimedec_limit_norm(nf, prime, BOUND);
    1353      231980 :     for (j = 1; j < lg(tabpr); j++)
    1354             :     {
    1355      115577 :       GEN pr  = gel(tabpr,j);
    1356      115577 :       if (condZ % p == 0 && idealval(nf, cond, pr))
    1357             :       {
    1358         462 :         gel(tmpray,j) = NULL; continue;
    1359             :       }
    1360      115115 :       vecsmalltrunc_append(R->L1, upowuu(p, pr_get_f(pr)));
    1361      115115 :       gel(tmpray,j) = gclone( isprincipalray(bnr, pr) );
    1362             :     }
    1363      116403 :     avma = av;
    1364      231980 :     for (k = 1; k < j; k++)
    1365             :     {
    1366      115577 :       if (!tmpray[k]) continue;
    1367      115115 :       vectrunc_append(R->L1ray, ZC_copy(gel(tmpray,k)));
    1368      115115 :       gunclone(gel(tmpray,k));
    1369             :     }
    1370             :   }
    1371         504 : }
    1372             : 
    1373             : static GEN /* cf polcoeff */
    1374      328440 : _sercoeff(GEN x, long n)
    1375             : {
    1376      328440 :   long i = n - valp(x);
    1377      328440 :   return (i < 0)? gen_0: gel(x,i+2);
    1378             : }
    1379             : 
    1380             : static void
    1381      328440 : affect_coeff(GEN q, long n, GEN y)
    1382             : {
    1383      328440 :   GEN x = _sercoeff(q,-n);
    1384      328440 :   if (x == gen_0) gel(y,n) = NULL; else affgr(x, gel(y,n));
    1385      328440 : }
    1386             : 
    1387             : typedef struct {
    1388             :   GEN c1, aij, bij, cS, cT, powracpi;
    1389             :   long i0, a,b,c, r, rc1, rc2;
    1390             : } ST_t;
    1391             : 
    1392             : /* compute the principal part at the integers s = 0, -1, -2, ..., -i0
    1393             :    of Gamma((s+1)/2)^a Gamma(s/2)^b Gamma(s)^c / (s - z) with z = 0 and 1 */
    1394             : /* NOTE: surely not the best way to do this, but it's fast enough! */
    1395             : static void
    1396         294 : ppgamma(ST_t *T, long prec)
    1397             : {
    1398             :   GEN eul, gam,gamun,gamdm, an,bn,cn_evn,cn_odd, x,x2,X,Y, cf, sqpi;
    1399             :   GEN p1, p2, aij, bij;
    1400         294 :   long a = T->a;
    1401         294 :   long b = T->b;
    1402         294 :   long c = T->c, r = T->r, i0 = T->i0;
    1403             :   long i,j, s,t;
    1404             :   pari_sp av;
    1405             : 
    1406         294 :   aij = cgetg(i0+1, t_VEC);
    1407         294 :   bij = cgetg(i0+1, t_VEC);
    1408       78120 :   for (i = 1; i <= i0; i++)
    1409             :   {
    1410       77826 :     gel(aij,i) = p1 = cgetg(r+1, t_VEC);
    1411       77826 :     gel(bij,i) = p2 = cgetg(r+1, t_VEC);
    1412       77826 :     for (j=1; j<=r; j++) { gel(p1,j) = cgetr(prec); gel(p2,j) = cgetr(prec); }
    1413             :   }
    1414         294 :   av = avma;
    1415             : 
    1416         294 :   x   = pol_x(0);
    1417         294 :   x2  = gmul2n(x, -1); /* x/2 */
    1418         294 :   eul = mpeuler(prec);
    1419         294 :   sqpi= sqrtr_abs(mppi(prec)); /* Gamma(1/2) */
    1420             : 
    1421             :   /* expansion of log(Gamma(u)) at u = 1 */
    1422         294 :   gamun = cgetg(r+3, t_SER);
    1423         294 :   gamun[1] = evalsigne(1) | _evalvalp(0) | evalvarn(0);
    1424         294 :   gel(gamun,2) = gen_0;
    1425         294 :   gel(gamun,3) = gneg(eul);
    1426         623 :   for (i = 2; i <= r; i++)
    1427         329 :     gel(gamun,i+2) = divrs(szeta(i,prec), odd(i)? -i: i);
    1428         294 :   gamun = gexp(gamun, prec); /* Gamma(1 + x) */
    1429         294 :   gam = gdiv(gamun,x); /* Gamma(x) */
    1430             : 
    1431             :   /* expansion of log(Gamma(u) / Gamma(1/2)) at u = 1/2 */
    1432         294 :   gamdm = cgetg(r+3, t_SER);
    1433         294 :   gamdm[1] = evalsigne(1) | _evalvalp(0) | evalvarn(0);
    1434         294 :   gel(gamdm,2) = gen_0;
    1435         294 :   gel(gamdm,3) = gneg(gadd(gmul2n(mplog2(prec), 1), eul));
    1436         623 :   for (i = 2; i <= r; i++)
    1437         329 :     gel(gamdm,i+2) = mulri(gel(gamun,i+2), subis(int2n(i), 1));
    1438         294 :   gamdm = gmul(sqpi, gexp(gamdm, prec)); /* Gamma(1/2 + x) */
    1439             : 
    1440             :  /* We simplify to get one of the following two expressions
    1441             :   * if (b > a) : sqrt{Pi}^a 2^{a-au} Gamma(u)^{a+c} Gamma(  u/2  )^{|b-a|}
    1442             :   * if (b <= a): sqrt{Pi}^b 2^{b-bu} Gamma(u)^{b+c} Gamma((u+1)/2)^{|b-a|} */
    1443         294 :   if (b > a)
    1444             :   {
    1445          49 :     t = a; s = b; X = x2; Y = gsub(x2,ghalf);
    1446          49 :     p1 = ser_unscale(gam, ghalf);
    1447          49 :     p2 = gdiv(ser_unscale(gamdm,ghalf), Y); /* Gamma((x-1)/2) */
    1448             :   }
    1449             :   else
    1450             :   {
    1451         245 :     t = b; s = a; X = gadd(x2,ghalf); Y = x2;
    1452         245 :     p1 = ser_unscale(gamdm,ghalf);
    1453         245 :     p2 = ser_unscale(gam,ghalf);
    1454             :   }
    1455         294 :   cf = powru(sqpi, t);
    1456         294 :   an = gpowgs(gpow(gen_2, gsubsg(1,x), prec), t); /* 2^{t-tx} */
    1457         294 :   bn = gpowgs(gam, t+c); /* Gamma(x)^{t+c} */
    1458         294 :   cn_evn = gpowgs(p1, s-t); /* Gamma(X)^{s-t} */
    1459         294 :   cn_odd = gpowgs(p2, s-t); /* Gamma(Y)^{s-t} */
    1460       39207 :   for (i = 0; i < i0/2; i++)
    1461             :   {
    1462       38913 :     GEN C1,q1, A1 = gel(aij,2*i+1), B1 = gel(bij,2*i+1);
    1463       38913 :     GEN C2,q2, A2 = gel(aij,2*i+2), B2 = gel(bij,2*i+2);
    1464             : 
    1465       38913 :     C1 = gmul(cf, gmul(bn, gmul(an, cn_evn)));
    1466       38913 :     p1 = gdiv(C1, gsubgs(x, 2*i));
    1467       38913 :     q1 = gdiv(C1, gsubgs(x, 2*i+1));
    1468             : 
    1469             :     /* an(x-u-1) = 2^t an(x-u) */
    1470       38913 :     an = gmul2n(an, t);
    1471             :     /* bn(x-u-1) = bn(x-u) / (x-u-1)^{t+c} */
    1472       38913 :     bn = gdiv(bn, gpowgs(gsubgs(x, 2*i+1), t+c));
    1473             : 
    1474       38913 :     C2 = gmul(cf, gmul(bn, gmul(an, cn_odd)));
    1475       38913 :     p2 = gdiv(C2, gsubgs(x, 2*i+1));
    1476       38913 :     q2 = gdiv(C2, gsubgs(x, 2*i+2));
    1477      121023 :     for (j = 1; j <= r; j++)
    1478             :     {
    1479       82110 :       affect_coeff(p1, j, A1); affect_coeff(q1, j, B1);
    1480       82110 :       affect_coeff(p2, j, A2); affect_coeff(q2, j, B2);
    1481             :     }
    1482             : 
    1483       38913 :     an = gmul2n(an, t);
    1484       38913 :     bn = gdiv(bn, gpowgs(gsubgs(x, 2*i+2), t+c));
    1485             :     /* cn_evn(x-2i-2) = cn_evn(x-2i)  / (X - (i+1))^{s-t} */
    1486             :     /* cn_odd(x-2i-3) = cn_odd(x-2i-1)/ (Y - (i+1))^{s-t} */
    1487       38913 :     cn_evn = gdiv(cn_evn, gpowgs(gsubgs(X,i+1), s-t));
    1488       38913 :     cn_odd = gdiv(cn_odd, gpowgs(gsubgs(Y,i+1), s-t));
    1489             :   }
    1490         294 :   T->aij = aij;
    1491         294 :   T->bij = bij; avma = av;
    1492         294 : }
    1493             : 
    1494             : static GEN
    1495        2324 : _cond(GEN dtcr) { return mkvec2(ch_cond(dtcr), ch_4(dtcr)); }
    1496             : 
    1497             : /* sort chars according to conductor */
    1498             : static GEN
    1499         700 : sortChars(GEN dataCR)
    1500             : {
    1501         700 :   const long cl = lg(dataCR) - 1;
    1502         700 :   GEN vCond  = cgetg(cl+1, t_VEC);
    1503         700 :   GEN CC     = cgetg(cl+1, t_VECSMALL);
    1504         700 :   GEN nvCond = cgetg(cl+1, t_VECSMALL);
    1505             :   long j,k, ncond;
    1506             :   GEN vChar;
    1507             : 
    1508         700 :   for (j = 1; j <= cl; j++) nvCond[j] = 0;
    1509             : 
    1510         700 :   ncond = 0;
    1511        3024 :   for (j = 1; j <= cl; j++)
    1512             :   {
    1513        2324 :     GEN cond = _cond(gel(dataCR,j));
    1514        2947 :     for (k = 1; k <= ncond; k++)
    1515        2058 :       if (gequal(cond, gel(vCond,k))) break;
    1516        2324 :     if (k > ncond) gel(vCond,++ncond) = cond;
    1517        2324 :     nvCond[k]++; CC[j] = k; /* char j has conductor number k */
    1518             :   }
    1519         700 :   vChar = cgetg(ncond+1, t_VEC);
    1520        1589 :   for (k = 1; k <= ncond; k++)
    1521             :   {
    1522         889 :     gel(vChar,k) = cgetg(nvCond[k]+1, t_VECSMALL);
    1523         889 :     nvCond[k] = 0;
    1524             :   }
    1525        3024 :   for (j = 1; j <= cl; j++)
    1526             :   {
    1527        2324 :     k = CC[j]; nvCond[k]++;
    1528        2324 :     mael(vChar,k,nvCond[k]) = j;
    1529             :   }
    1530         700 :   return vChar;
    1531             : }
    1532             : 
    1533             : /* Given W(chi), S(chi) and T(chi), return L(1, chi) if fl & 1, else
    1534             :    [r(chi), c(chi)] where L(s, chi) ~ c(chi) s^r(chi) at s = 0.
    1535             :    If fl & 2, adjust the value to get L_S(s, chi). */
    1536             : static GEN
    1537        1337 : GetValue(GEN dtcr, GEN W, GEN S, GEN T, long fl, long prec)
    1538             : {
    1539        1337 :   pari_sp av = avma;
    1540             :   GEN cf, z, p1;
    1541             :   long q, b, c, r;
    1542        1337 :   int isreal = (chi_get_deg(ch_CHI0(dtcr)) <= 2);
    1543             : 
    1544        1337 :   p1 = ch_4(dtcr);
    1545        1337 :   q = p1[1];
    1546        1337 :   b = p1[2];
    1547        1337 :   c = p1[3];
    1548             : 
    1549        1337 :   if (fl & 1)
    1550             :   { /* S(chi) + W(chi).T(chi)) / (C(chi) sqrt(Pi)^{r1 - q}) */
    1551         196 :     cf = gmul(ch_C(dtcr), powruhalf(mppi(prec), b));
    1552             : 
    1553         196 :     z = gadd(S, gmul(W, T));
    1554         196 :     if (isreal) z = real_i(z);
    1555         196 :     z = gdiv(z, cf);
    1556         196 :     if (fl & 2) z = gmul(z, ComputeAChi(dtcr, &r, 1, prec));
    1557             :   }
    1558             :   else
    1559             :   { /* (W(chi).S(conj(chi)) + T(chi)) / (sqrt(Pi)^q 2^{r1 - q}) */
    1560        1141 :     cf = gmul2n(powruhalf(mppi(prec), q), b);
    1561             : 
    1562        1141 :     z = gadd(gmul(W, gconj(S)), gconj(T));
    1563        1141 :     if (isreal) z = real_i(z);
    1564        1141 :     z = gdiv(z, cf); r = 0;
    1565        1141 :     if (fl & 2) z = gmul(z, ComputeAChi(dtcr, &r, 0, prec));
    1566        1141 :     z = mkvec2(utoi(b + c + r), z);
    1567             :   }
    1568        1337 :   return gerepilecopy(av, z);
    1569             : }
    1570             : 
    1571             : /* return the order and the first non-zero term of L(s, chi0)
    1572             :    at s = 0. If flag != 0, adjust the value to get L_S(s, chi0). */
    1573             : static GEN
    1574          21 : GetValue1(GEN bnr, long flag, long prec)
    1575             : {
    1576          21 :   GEN bnf = checkbnf(bnr), nf = bnf_get_nf(bnf);
    1577             :   GEN h, R, c, diff;
    1578             :   long i, l, r, r1, r2;
    1579          21 :   pari_sp av = avma;
    1580             : 
    1581          21 :   nf_get_sign(nf, &r1,&r2);
    1582          21 :   h = bnf_get_no(bnf);
    1583          21 :   R = bnf_get_reg(bnf);
    1584             : 
    1585          21 :   c = gneg_i(gdivgs(mpmul(h, R), bnf_get_tuN(bnf)));
    1586          21 :   r = r1 + r2 - 1;
    1587             : 
    1588          21 :   if (flag)
    1589             :   {
    1590           0 :     diff = divcond(bnr);
    1591           0 :     l = lg(diff) - 1; r += l;
    1592           0 :     for (i = 1; i <= l; i++)
    1593           0 :       c = gmul(c, glog(pr_norm(gel(diff,i)), prec));
    1594             :   }
    1595          21 :   return gerepilecopy(av, mkvec2(stoi(r), c));
    1596             : }
    1597             : 
    1598             : /********************************************************************/
    1599             : /*                6th part: recover the coefficients                */
    1600             : /********************************************************************/
    1601             : static long
    1602        2505 : TestOne(GEN plg, RC_data *d)
    1603             : {
    1604        2505 :   long j, v = d->v;
    1605        2505 :   GEN z = gsub(d->beta, gel(plg,v));
    1606        2505 :   if (expo(z) >= d->G) return 0;
    1607        6807 :   for (j = 1; j < lg(plg); j++)
    1608        4750 :     if (j != v && mpcmp(d->B, mpabs(gel(plg,j))) < 0) return 0;
    1609        2057 :   return 1;
    1610             : }
    1611             : 
    1612             : static GEN
    1613         503 : chk_reccoeff_init(FP_chk_fun *chk, GEN r, GEN mat)
    1614             : {
    1615         503 :   RC_data *d = (RC_data*)chk->data;
    1616         503 :   (void)r; d->U = mat; return d->nB;
    1617             : }
    1618             : 
    1619             : static GEN
    1620         530 : chk_reccoeff(void *data, GEN x)
    1621             : {
    1622         530 :   RC_data *d = (RC_data*)data;
    1623         530 :   GEN v = gmul(d->U, x), z = gel(v,1);
    1624             : 
    1625         530 :   if (!gequal1(z)) return NULL;
    1626         530 :   *++v = evaltyp(t_COL) | evallg( lg(d->M) );
    1627         530 :   if (TestOne(gmul(d->M, v), d)) return v;
    1628           0 :   return NULL;
    1629             : }
    1630             : 
    1631             : /* Using Cohen's method */
    1632             : static GEN
    1633         503 : RecCoeff3(GEN nf, RC_data *d, long prec)
    1634             : {
    1635             :   GEN A, M, nB, cand, p1, B2, C2, tB, beta2, nf2, Bd;
    1636         503 :   GEN beta = d->beta, B = d->B;
    1637         503 :   long N = d->N, v = d->v, e, BIG;
    1638         503 :   long i, j, k, ct = 0, prec2;
    1639         503 :   FP_chk_fun chk = { &chk_reccoeff, &chk_reccoeff_init, NULL, NULL, 0 };
    1640         503 :   chk.data = (void*)d;
    1641             : 
    1642         503 :   d->G = minss(-10, -prec2nbits(prec) >> 4);
    1643         503 :   BIG = maxss(32, -2*d->G);
    1644         503 :   tB  = sqrtnr(real2n(BIG-N,DEFAULTPREC), N-1);
    1645         503 :   Bd  = grndtoi(gmin(B, tB), &e);
    1646         503 :   if (e > 0) return NULL; /* failure */
    1647         503 :   Bd = addis(Bd, 1);
    1648         503 :   prec2 = nbits2prec( expi(Bd) + 192 );
    1649         503 :   prec2 = maxss(precdbl(prec), prec2);
    1650         503 :   B2 = sqri(Bd);
    1651         503 :   C2 = shifti(B2, BIG<<1);
    1652             : 
    1653         503 : LABrcf: ct++;
    1654         503 :   beta2 = gprec_w(beta, prec2);
    1655         503 :   nf2 = nfnewprec_shallow(nf, prec2);
    1656         503 :   d->M = M = nf_get_M(nf2);
    1657             : 
    1658         503 :   A = cgetg(N+2, t_MAT);
    1659         503 :   for (i = 1; i <= N+1; i++) gel(A,i) = cgetg(N+2, t_COL);
    1660             : 
    1661         503 :   gcoeff(A, 1, 1) = gadd(gmul(C2, gsqr(beta2)), B2);
    1662        1607 :   for (j = 2; j <= N+1; j++)
    1663             :   {
    1664        1104 :     p1 = gmul(C2, gmul(gneg_i(beta2), gcoeff(M, v, j-1)));
    1665        1104 :     gcoeff(A, 1, j) = gcoeff(A, j, 1) = p1;
    1666             :   }
    1667        1607 :   for (i = 2; i <= N+1; i++)
    1668        2907 :     for (j = i; j <= N+1; j++)
    1669             :     {
    1670        1803 :       p1 = gen_0;
    1671        5997 :       for (k = 1; k <= N; k++)
    1672             :       {
    1673        4194 :         GEN p2 = gmul(gcoeff(M, k, j-1), gcoeff(M, k, i-1));
    1674        4194 :         if (k == v) p2 = gmul(C2, p2);
    1675        4194 :         p1 = gadd(p1,p2);
    1676             :       }
    1677        1803 :       gcoeff(A, i, j) = gcoeff(A, j, i) = p1;
    1678             :     }
    1679             : 
    1680         503 :   nB = mului(N+1, B2);
    1681         503 :   d->nB = nB;
    1682         503 :   cand = fincke_pohst(A, nB, -1, prec2, &chk);
    1683             : 
    1684         503 :   if (!cand)
    1685             :   {
    1686           0 :     if (ct > 3) return NULL;
    1687           0 :     prec2 = precdbl(prec2);
    1688           0 :     if (DEBUGLEVEL>1) pari_warn(warnprec,"RecCoeff", prec2);
    1689           0 :     goto LABrcf;
    1690             :   }
    1691             : 
    1692         503 :   cand = gel(cand,1);
    1693         503 :   if (lg(cand) == 2) return gel(cand,1);
    1694             : 
    1695         238 :   if (DEBUGLEVEL>1) err_printf("RecCoeff3: no solution found!\n");
    1696         238 :   return NULL;
    1697             : }
    1698             : 
    1699             : /* Using linear dependance relations */
    1700             : static GEN
    1701        2030 : RecCoeff2(GEN nf,  RC_data *d,  long prec)
    1702             : {
    1703             :   pari_sp av;
    1704        2030 :   GEN vec, M = nf_get_M(nf), beta = d->beta;
    1705        2030 :   long i, imin, imax, lM = lg(M);
    1706             : 
    1707        2030 :   d->G = minss(-20, -prec2nbits(prec) >> 4);
    1708             : 
    1709        2030 :   vec  = shallowconcat(mkvec(gneg(beta)), row(M, d->v));
    1710        2030 :   imin = (long)prec2nbits_mul(prec, .225);
    1711        2030 :   imax = (long)prec2nbits_mul(prec, .315);
    1712             : 
    1713        2030 :   av = avma;
    1714        2567 :   for (i = imax; i >= imin; i-=16, avma = av)
    1715             :   {
    1716             :     long e;
    1717        2064 :     GEN v = lindep2(vec, i), z = gel(v,1);
    1718        2064 :     if (!signe(z)) continue;
    1719        1975 :     *++v = evaltyp(t_COL) | evallg(lM);
    1720        1975 :     v = grndtoi(gdiv(v, z), &e);
    1721        1975 :     if (e > 0) break;
    1722        1975 :     if (TestOne(gmul(M, v), d)) return v;
    1723             :   }
    1724             :   /* failure */
    1725         503 :   return RecCoeff3(nf,d,prec);
    1726             : }
    1727             : 
    1728             : /* Attempts to find a polynomial with coefficients in nf such that
    1729             :    its coefficients are close to those of pol at the place v and
    1730             :    less than B at all the other places */
    1731             : static GEN
    1732         581 : RecCoeff(GEN nf,  GEN pol,  long v, long prec)
    1733             : {
    1734         581 :   long j, md, cl = degpol(pol);
    1735         581 :   pari_sp av = avma;
    1736             :   RC_data d;
    1737             : 
    1738             :   /* if precision(pol) is too low, abort */
    1739        3717 :   for (j = 2; j <= cl+1; j++)
    1740             :   {
    1741        3157 :     GEN t = gel(pol, j);
    1742        3157 :     if (prec2nbits(gprecision(t)) - gexpo(t) < 34) return NULL;
    1743             :   }
    1744             : 
    1745         560 :   md = cl/2;
    1746         560 :   pol = leafcopy(pol);
    1747             : 
    1748         560 :   d.N = nf_get_degree(nf);
    1749         560 :   d.v = v;
    1750             : 
    1751        2352 :   for (j = 1; j <= cl; j++)
    1752             :   { /* start with the coefficients in the middle,
    1753             :        since they are the harder to recognize! */
    1754        2030 :     long cf = md + (j%2? j/2: -j/2);
    1755        2030 :     GEN t, bound = shifti(binomial(utoipos(cl), cf), cl-cf);
    1756             : 
    1757        2030 :     if (DEBUGLEVEL>1) err_printf("RecCoeff (cf = %ld, B = %Ps)\n", cf, bound);
    1758        2030 :     d.beta = real_i( gel(pol,cf+2) );
    1759        2030 :     d.B    = bound;
    1760        2030 :     if (! (t = RecCoeff2(nf, &d, prec)) ) return NULL;
    1761        1792 :     gel(pol, cf+2) = coltoalg(nf,t);
    1762             :   }
    1763         322 :   gel(pol,cl+2) = gen_1;
    1764         322 :   return gerepilecopy(av, pol);
    1765             : }
    1766             : 
    1767             : /* an[q * i] *= chi for all (i,p)=1 */
    1768             : static void
    1769      137187 : an_mul(int **an, long p, long q, long n, long deg, GEN chi, int **reduc)
    1770             : {
    1771             :   pari_sp av;
    1772             :   long c,i;
    1773             :   int *T;
    1774             : 
    1775      274374 :   if (gequal1(chi)) return;
    1776      126469 :   av = avma;
    1777      126469 :   T = (int*)new_chunk(deg); Polmod2Coeff(T,chi, deg);
    1778     2463428 :   for (c = 1, i = q; i <= n; i += q, c++)
    1779     2336959 :     if (c == p) c = 0; else MulCoeff(an[i], T, reduc, deg);
    1780      126469 :   avma = av;
    1781             : }
    1782             : /* an[q * i] = 0 for all (i,p)=1 */
    1783             : static void
    1784      125225 : an_set0_coprime(int **an, long p, long q, long n, long deg)
    1785             : {
    1786             :   long c,i;
    1787     1600719 :   for (c = 1, i = q; i <= n; i += q, c++)
    1788     1475494 :     if (c == p) c = 0; else _0toCoeff(an[i], deg);
    1789      125225 : }
    1790             : /* an[q * i] = 0 for all i */
    1791             : static void
    1792         280 : an_set0(int **an, long p, long n, long deg)
    1793             : {
    1794             :   long i;
    1795         280 :   for (i = p; i <= n; i += p) _0toCoeff(an[i], deg);
    1796         280 : }
    1797             : 
    1798             : /* compute the coefficients an for the quadratic case */
    1799             : static int**
    1800         791 : computean(GEN dtcr, LISTray *R, long n, long deg)
    1801             : {
    1802         791 :   pari_sp av = avma, av2;
    1803             :   long i, p, q, condZ, l;
    1804             :   int **an, **reduc;
    1805             :   GEN L, chi, chi1;
    1806             :   CHI_t C;
    1807             : 
    1808         791 :   init_CHI_alg(&C, ch_CHI(dtcr));
    1809         791 :   condZ= R->condZ;
    1810             : 
    1811         791 :   an = InitMatAn(n, deg, 1);
    1812         791 :   reduc = InitReduction(C.ord, deg);
    1813         791 :   av2 = avma;
    1814             : 
    1815             :   /* all pr | p divide cond */
    1816         791 :   L = R->L0; l = lg(L);
    1817         791 :   for (i=1; i<l; i++) an_set0(an,L[i],n,deg);
    1818             : 
    1819             :   /* 1 prime of degree 2 */
    1820         791 :   L = R->L2; l = lg(L);
    1821      124341 :   for (i=1; i<l; i++, avma = av2)
    1822             :   {
    1823      123550 :     p = L[i];
    1824      123550 :     if (condZ == 1) chi = C.val[0]; /* 1 */
    1825      123396 :     else            chi = CHI_eval(&C, gel(R->rayZ, p%condZ));
    1826      123550 :     chi1 = chi;
    1827      123550 :     for (q=p;;)
    1828             :     {
    1829      125225 :       an_set0_coprime(an, p,q,n,deg); /* v_p(q) odd */
    1830      125225 :       if (! (q = next_pow(q,p, n)) ) break;
    1831             : 
    1832        5376 :       an_mul(an,p,q,n,deg,chi,reduc);
    1833        5376 :       if (! (q = next_pow(q,p, n)) ) break;
    1834        1675 :       chi = gmul(chi, chi1);
    1835        1675 :     }
    1836             :   }
    1837             : 
    1838             :   /* 1 prime of degree 1 */
    1839         791 :   L = R->L1; l = lg(L);
    1840        2394 :   for (i=1; i<l; i++, avma = av2)
    1841             :   {
    1842        1603 :     p = L[i];
    1843        1603 :     chi = CHI_eval(&C, gel(R->L1ray,i));
    1844        1603 :     chi1 = chi;
    1845        1603 :     for(q=p;;)
    1846             :     {
    1847        7693 :       an_mul(an,p,q,n,deg,chi,reduc);
    1848        7693 :       if (! (q = next_pow(q,p, n)) ) break;
    1849        6090 :       chi = gmul(chi, chi1);
    1850        6090 :     }
    1851             :   }
    1852             : 
    1853             :   /* 2 primes of degree 1 */
    1854         791 :   L = R->L11; l = lg(L);
    1855      117489 :   for (i=1; i<l; i++, avma = av2)
    1856             :   {
    1857             :     GEN ray1, ray2, chi11, chi12, chi2;
    1858             : 
    1859      116698 :     p = L[i]; ray1 = gel(R->L11ray,i); /* use pr1 pr2 = (p) */
    1860      116698 :     if (condZ == 1)
    1861         112 :       ray2 = ZC_neg(ray1);
    1862             :     else
    1863      116586 :       ray2 = ZC_sub(gel(R->rayZ, p%condZ),  ray1);
    1864      116698 :     chi11 = CHI_eval(&C, ray1);
    1865      116698 :     chi12 = CHI_eval(&C, ray2);
    1866             : 
    1867      116698 :     chi1 = gadd(chi11, chi12);
    1868      116698 :     chi2 = chi12;
    1869      116698 :     for(q=p;;)
    1870             :     {
    1871      124118 :       an_mul(an,p,q,n,deg,chi1,reduc);
    1872      124118 :       if (! (q = next_pow(q,p, n)) ) break;
    1873        7420 :       chi2 = gmul(chi2, chi12);
    1874        7420 :       chi1 = gadd(chi2, gmul(chi1, chi11));
    1875        7420 :     }
    1876             :   }
    1877             : 
    1878         791 :   CorrectCoeff(dtcr, an, reduc, n, deg);
    1879         791 :   FreeMat(reduc, deg-1);
    1880         791 :   avma = av; return an;
    1881             : }
    1882             : 
    1883             : /* return the vector of A^i/i for i = 1...n */
    1884             : static GEN
    1885         252 : mpvecpowdiv(GEN A, long n)
    1886             : {
    1887         252 :   pari_sp av = avma;
    1888             :   long i;
    1889         252 :   GEN v = powersr(A, n);
    1890         252 :   GEN w = cgetg(n+1, t_VEC);
    1891         252 :   gel(w,1) = rcopy(gel(v,2));
    1892         252 :   for (i=2; i<=n; i++) gel(w,i) = divru(gel(v,i+1), i);
    1893         252 :   return gerepileupto(av, w);
    1894             : }
    1895             : 
    1896             : static void GetST0(GEN bnr, GEN *pS, GEN *pT, GEN dataCR, GEN vChar, long prec);
    1897             : 
    1898             : /* compute S and T for the quadratic case. The following cases (cs) are:
    1899             :    1) bnr complex;
    1900             :    2) bnr real and no infinite place divide cond_chi (TBD);
    1901             :    3) bnr real and one infinite place divide cond_chi;
    1902             :    4) bnr real and both infinite places divide cond_chi (TBD) */
    1903             : static void
    1904         252 : QuadGetST(GEN bnr, GEN *pS, GEN *pT, GEN dataCR, GEN vChar, long prec)
    1905             : {
    1906         252 :   pari_sp av = avma, av1, av2;
    1907             :   long ncond, n, j, k, n0;
    1908         252 :   GEN N0, C, T = *pT, S = *pS, an, degs, cs;
    1909             :   LISTray LIST;
    1910             : 
    1911             :   /* initializations */
    1912         252 :   degs = GetDeg(dataCR);
    1913         252 :   ncond = lg(vChar)-1;
    1914         252 :   C    = cgetg(ncond+1, t_VEC);
    1915         252 :   N0   = cgetg(ncond+1, t_VECSMALL);
    1916         252 :   cs   = cgetg(ncond+1, t_VECSMALL);
    1917         252 :   n0 = 0;
    1918         525 :   for (j = 1; j <= ncond; j++)
    1919             :   {
    1920             :     /* FIXME: make sure that this value of c is correct for the general case */
    1921             :     long r1, r2, q;
    1922         308 :     GEN dtcr = gel(dataCR, mael(vChar,j,1)), p1 = ch_4(dtcr), c = ch_C(dtcr);
    1923             : 
    1924         308 :     gel(C,j) = c;
    1925         308 :     q = p1[1];
    1926             : 
    1927         308 :     nf_get_sign(bnr_get_nf(ch_bnr(dtcr)), &r1, &r2);
    1928         308 :     if (r1 == 2) /* real quadratic */
    1929             :     {
    1930         294 :       cs[j] = 2 + q;
    1931             :       /* FIXME:
    1932             :          make sure that this value of N0 is correct for the general case */
    1933         294 :       N0[j] = (long)prec2nbits_mul(prec, 0.35 * gtodouble(c));
    1934         294 :       if (cs[j] == 2 || cs[j] == 4) /* NOT IMPLEMENTED YET */
    1935             :       {
    1936          35 :         GetST0(bnr, pS, pT, dataCR, vChar, prec);
    1937         287 :         return;
    1938             :       }
    1939             :     }
    1940             :     else /* complex quadratic */
    1941             :     {
    1942          14 :       cs[j] = 1;
    1943          14 :       N0[j] = (long)prec2nbits_mul(prec, 0.7 * gtodouble(c));
    1944             :     }
    1945         273 :     if (n0 < N0[j]) n0 = N0[j];
    1946             :   }
    1947         217 :   if (DEBUGLEVEL>1) err_printf("N0 = %ld\n", n0);
    1948         217 :   InitPrimesQuad(bnr, n0, &LIST);
    1949             : 
    1950         217 :   av1 = avma;
    1951             :   /* loop over conductors */
    1952         469 :   for (j = 1; j <= ncond; j++)
    1953             :   {
    1954         252 :     GEN c0 = gel(C,j), c1 = divur(1, c0), c2 = divur(2, c0);
    1955         252 :     GEN ec1 = mpexp(c1), ec2 = mpexp(c2), LChar = gel(vChar,j);
    1956             :     GEN vf0, vf1, cf0, cf1;
    1957         252 :     const long nChar = lg(LChar)-1, NN = N0[j];
    1958             : 
    1959         252 :     if (DEBUGLEVEL>1)
    1960           0 :       err_printf("* conductor no %ld/%ld (N = %ld)\n\tInit: ", j,ncond,NN);
    1961         252 :     if (realprec(ec1) > prec) ec1 = rtor(ec1, prec);
    1962         252 :     if (realprec(ec2) > prec) ec2 = rtor(ec2, prec);
    1963         252 :     switch(cs[j])
    1964             :     {
    1965             :     case 1:
    1966          14 :       cf0 = gen_1;
    1967          14 :       cf1 = c0;
    1968          14 :       vf0 = mpveceint1(rtor(c1, prec), ec1, NN);
    1969          14 :       vf1 = mpvecpowdiv(invr(ec1), NN); break;
    1970             : 
    1971             :     case 3:
    1972         238 :       cf0 = sqrtr(mppi(prec));
    1973         238 :       cf1 = gmul2n(cf0, 1);
    1974         238 :       cf0 = gmul(cf0, c0);
    1975         238 :       vf0 = mpvecpowdiv(invr(ec2), NN);
    1976         238 :       vf1 = mpveceint1(rtor(c2, prec), ec2, NN); break;
    1977             : 
    1978             :     default:
    1979           0 :       cf0 = cf1 = NULL; /* FIXME: not implemented */
    1980           0 :       vf0 = vf1 = NULL;
    1981             :     }
    1982        1050 :     for (k = 1; k <= nChar; k++)
    1983             :     {
    1984         798 :       const long t = LChar[k], d = degs[t];
    1985         798 :       const GEN dtcr = gel(dataCR, t), z = gel(ch_CHI(dtcr), 2);
    1986         798 :       GEN p1 = gen_0, p2 = gen_0;
    1987             :       int **matan;
    1988         798 :       long c = 0;
    1989             : 
    1990         798 :       if (DEBUGLEVEL>1)
    1991           0 :         err_printf("\tcharacter no: %ld (%ld/%ld)\n", t,k,nChar);
    1992         798 :       if (isintzero( ch_comp(gel(dataCR, t)) ))
    1993             :       {
    1994           7 :         if (DEBUGLEVEL>1) err_printf("\t  no need to compute this character\n");
    1995           7 :         continue;
    1996             :       }
    1997         791 :       av2 = avma;
    1998         791 :       matan = computean(gel(dataCR,t), &LIST, NN, d);
    1999     1599884 :       for (n = 1; n <= NN; n++)
    2000     1599093 :         if ((an = EvalCoeff(z, matan[n], d)))
    2001             :         {
    2002      334501 :           p1 = gadd(p1, gmul(an, gel(vf0,n)));
    2003      334501 :           p2 = gadd(p2, gmul(an, gel(vf1,n)));
    2004      334501 :           if (++c == 256) { gerepileall(av2,2, &p1,&p2); c = 0; }
    2005             :         }
    2006         791 :       gaffect(gmul(cf0, p1), gel(S,t));
    2007         791 :       gaffect(gmul(cf1,  gconj(p2)), gel(T,t));
    2008         791 :       FreeMat(matan,NN); avma = av2;
    2009             :     }
    2010         252 :     if (DEBUGLEVEL>1) err_printf("\n");
    2011         252 :     avma = av1;
    2012             :   }
    2013         217 :   avma = av;
    2014             : }
    2015             : 
    2016             : /* s += t*u. All 3 of them t_REAL, except we allow s or u = NULL (for 0) */
    2017             : static GEN
    2018    49382186 : _addmulrr(GEN s, GEN t, GEN u)
    2019             : {
    2020    49382186 :   if (u)
    2021             :   {
    2022    49119427 :     GEN v = mulrr(t, u);
    2023    49119427 :     return s? addrr(s, v): v;
    2024             :   }
    2025      262759 :   return s;
    2026             : }
    2027             : /* s += t. Both real, except we allow s or t = NULL (for exact 0) */
    2028             : static GEN
    2029    99692166 : _addrr(GEN s, GEN t)
    2030    99692166 : { return t? (s? addrr(s, t): t) : s; }
    2031             : 
    2032             : /* S & T for the general case. This is time-critical: optimize */
    2033             : static void
    2034      503195 : get_cS_cT(ST_t *T, long n)
    2035             : {
    2036             :   pari_sp av;
    2037             :   GEN csurn, nsurc, lncsurn, A, B, s, t, Z, aij, bij;
    2038             :   long i, j, r, i0;
    2039             : 
    2040     1006390 :   if (T->cS[n]) return;
    2041             : 
    2042      239106 :   av = avma;
    2043      239106 :   aij = T->aij; i0= T->i0;
    2044      239106 :   bij = T->bij; r = T->r;
    2045      239106 :   Z = cgetg(r+1, t_VEC);
    2046      239106 :   gel(Z,1) = NULL; /* unused */
    2047             : 
    2048      239106 :   csurn = divru(T->c1, n);
    2049      239106 :   nsurc = invr(csurn);
    2050      239106 :   lncsurn = logr_abs(csurn);
    2051             : 
    2052      239106 :   if (r > 1)
    2053             :   {
    2054      239043 :     gel(Z,2) = lncsurn; /* r >= 2 */
    2055      243138 :     for (i = 3; i <= r; i++)
    2056        4095 :       gel(Z,i) = divru(mulrr(gel(Z,i-1), lncsurn), i-1);
    2057             :     /* Z[i] = ln^(i-1)(c1/n) / (i-1)! */
    2058             :   }
    2059             : 
    2060             :   /* i = i0 */
    2061      239106 :     A = gel(aij,i0); t = _addrr(NULL, gel(A,1));
    2062      239106 :     B = gel(bij,i0); s = _addrr(NULL, gel(B,1));
    2063      482244 :     for (j = 2; j <= r; j++)
    2064             :     {
    2065      243138 :       s = _addmulrr(s, gel(Z,j),gel(B,j));
    2066      243138 :       t = _addmulrr(t, gel(Z,j),gel(A,j));
    2067             :     }
    2068    49487424 :   for (i = i0 - 1; i > 1; i--)
    2069             :   {
    2070    49248318 :     A = gel(aij,i); if (t) t = mulrr(t, nsurc);
    2071    49248318 :     B = gel(bij,i); if (s) s = mulrr(s, nsurc);
    2072    73453135 :     for (j = odd(i)? T->rc2: T->rc1; j > 1; j--)
    2073             :     {
    2074    24204817 :       s = _addmulrr(s, gel(Z,j),gel(B,j));
    2075    24204817 :       t = _addmulrr(t, gel(Z,j),gel(A,j));
    2076             :     }
    2077    49248318 :     s = _addrr(s, gel(B,1));
    2078    49248318 :     t = _addrr(t, gel(A,1));
    2079             :   }
    2080             :   /* i = 1 */
    2081      239106 :     A = gel(aij,1); if (t) t = mulrr(t, nsurc);
    2082      239106 :     B = gel(bij,1); if (s) s = mulrr(s, nsurc);
    2083      239106 :     s = _addrr(s, gel(B,1));
    2084      239106 :     t = _addrr(t, gel(A,1));
    2085      482244 :     for (j = 2; j <= r; j++)
    2086             :     {
    2087      243138 :       s = _addmulrr(s, gel(Z,j),gel(B,j));
    2088      243138 :       t = _addmulrr(t, gel(Z,j),gel(A,j));
    2089             :     }
    2090      239106 :   s = _addrr(s, T->b? mulrr(csurn, gel(T->powracpi,T->b+1)): csurn);
    2091      239106 :   if (!s) s = gen_0;
    2092      239106 :   if (!t) t = gen_0;
    2093      239106 :   gel(T->cS,n) = gclone(s);
    2094      239106 :   gel(T->cT,n) = gclone(t); avma = av;
    2095             : }
    2096             : 
    2097             : static void
    2098         476 : clear_cScT(ST_t *T, long N)
    2099             : {
    2100         476 :   GEN cS = T->cS, cT = T->cT;
    2101             :   long i;
    2102     1517068 :   for (i=1; i<=N; i++)
    2103     1516592 :     if (cS[i]) {
    2104      239106 :       gunclone(gel(cS,i));
    2105      239106 :       gunclone(gel(cT,i)); gel(cS,i) = gel(cT,i) = NULL;
    2106             :     }
    2107         476 : }
    2108             : 
    2109             : static void
    2110         294 : init_cScT(ST_t *T, GEN dtcr, long N, long prec)
    2111             : {
    2112         294 :   GEN p1 = ch_4(dtcr);
    2113         294 :   T->a = p1[1];
    2114         294 :   T->b = p1[2];
    2115         294 :   T->c = p1[3];
    2116         294 :   T->rc1 = T->a + T->c;
    2117         294 :   T->rc2 = T->b + T->c;
    2118         294 :   T->r   = maxss(T->rc2+1, T->rc1); /* >= 2 */
    2119         294 :   ppgamma(T, prec);
    2120         294 :   clear_cScT(T, N);
    2121         294 : }
    2122             : 
    2123             : /* return a t_REAL */
    2124             : static GEN
    2125         504 : zeta_get_limx(long r1, long r2, long bit)
    2126             : {
    2127         504 :   pari_sp av = avma;
    2128             :   GEN p1, p2, c0, c1, A0;
    2129         504 :   long r = r1 + r2, N = r + r2;
    2130             : 
    2131             :   /* c1 = N 2^(-2r2 / N) */
    2132         504 :   c1 = mulrs(powrfrac(real2n(1, DEFAULTPREC), -2*r2, N), N);
    2133             : 
    2134         504 :   p1 = powru(Pi2n(1, DEFAULTPREC), r - 1);
    2135         504 :   p2 = mulir(powuu(N,r), p1); shiftr_inplace(p2, -r2);
    2136         504 :   c0 = sqrtr( divrr(p2, powru(c1, r+1)) );
    2137             : 
    2138         504 :   A0 = logr_abs( gmul2n(c0, bit) ); p2 = divrr(A0, c1);
    2139         504 :   p1 = divrr(mulur(N*(r+1), logr_abs(p2)), addsr(2*(r+1), gmul2n(A0,2)));
    2140         504 :   return gerepileuptoleaf(av, divrr(addrs(p1, 1), powruhalf(p2, N)));
    2141             : }
    2142             : /* N_0 = floor( C_K / limx ). Large */
    2143             : static long
    2144         616 : zeta_get_N0(GEN C,  GEN limx)
    2145             : {
    2146             :   long e;
    2147         616 :   pari_sp av = avma;
    2148         616 :   GEN z = gcvtoi(gdiv(C, limx), &e); /* avoid truncation error */
    2149         616 :   if (e >= 0 || is_bigint(z))
    2150           0 :     pari_err_OVERFLOW("zeta_get_N0 [need too many primes]");
    2151         616 :   if (DEBUGLEVEL>1) err_printf("\ninitzeta: N0 = %Ps\n", z);
    2152         616 :   avma = av; return itos(z);
    2153             : }
    2154             : 
    2155             : /* even i such that limx^i ( (i\2)! )^r1 ( i! )^r2 ~ B. Small. */
    2156             : static long
    2157         182 : get_i0(long r1, long r2, GEN B, GEN limx)
    2158             : {
    2159         182 :   long imin = 1, imax = 1400;
    2160        2002 :   while(imax - imin >= 4)
    2161             :   {
    2162        1638 :     long i = (imax + imin) >> 1;
    2163        1638 :     GEN t = powru(limx, i);
    2164        1638 :     if (!r1)      t = mulrr(t, powru(mpfactr(i  , DEFAULTPREC), r2));
    2165        1638 :     else if (!r2) t = mulrr(t, powru(mpfactr(i/2, DEFAULTPREC), r1));
    2166             :     else {
    2167           0 :       GEN u1 = mpfactr(i/2, DEFAULTPREC);
    2168           0 :       GEN u2 = mpfactr(i,   DEFAULTPREC);
    2169           0 :       if (r1 == r2) t = mulrr(t, powru(mulrr(u1,u2), r1));
    2170           0 :       else t = mulrr(t, mulrr(powru(u1,r1), powru(u2,r2)));
    2171             :     }
    2172        1638 :     if (mpcmp(t, B) >= 0) imax = i; else imin = i;
    2173             :   }
    2174         182 :   return imax & ~1; /* make it even */
    2175             : }
    2176             : /* assume limx = zeta_get_limx(r1, r2, bit), a t_REAL */
    2177             : static long
    2178         182 : zeta_get_i0(long r1, long r2, long bit, GEN limx)
    2179             : {
    2180         182 :   pari_sp av = avma;
    2181         182 :   GEN B = gmul(sqrtr( divrr(powrs(mppi(DEFAULTPREC), r2-3), limx) ),
    2182             :                gmul2n(powuu(5, r1), bit + r2));
    2183         182 :   long i0 = get_i0(r1, r2, B, limx);
    2184         182 :   if (DEBUGLEVEL>1) { err_printf("i0 = %ld\n",i0); err_flush(); }
    2185         182 :   avma = av; return i0;
    2186             : }
    2187             : 
    2188             : static void
    2189         182 : GetST0(GEN bnr, GEN *pS, GEN *pT, GEN dataCR, GEN vChar, long prec)
    2190             : {
    2191         182 :   pari_sp av = avma, av1, av2;
    2192             :   long ncond, n, j, k, jc, n0, prec2, i0, r1, r2;
    2193         182 :   GEN nf = checknf(bnr), T = *pT, S = *pS;
    2194             :   GEN N0, C, an, degs, limx;
    2195             :   LISTray LIST;
    2196             :   ST_t cScT;
    2197             : 
    2198             :   /* initializations */
    2199         182 :   degs = GetDeg(dataCR);
    2200         182 :   ncond = lg(vChar)-1;
    2201         182 :   nf_get_sign(nf,&r1,&r2);
    2202             : 
    2203         182 :   C  = cgetg(ncond+1, t_VEC);
    2204         182 :   N0 = cgetg(ncond+1, t_VECSMALL);
    2205         182 :   n0 = 0;
    2206         182 :   limx = zeta_get_limx(r1, r2, prec2nbits(prec));
    2207         476 :   for (j = 1; j <= ncond; j++)
    2208             :   {
    2209         294 :     GEN dtcr = gel(dataCR, mael(vChar,j,1)), c = ch_C(dtcr);
    2210         294 :     gel(C,j) = c;
    2211         294 :     N0[j] = zeta_get_N0(c, limx);
    2212         294 :     if (n0 < N0[j]) n0  = N0[j];
    2213             :   }
    2214         182 :   i0 = zeta_get_i0(r1, r2, prec2nbits(prec), limx);
    2215         182 :   InitPrimes(bnr, n0, &LIST);
    2216             : 
    2217         182 :   prec2 = precdbl(prec) + EXTRA_PREC;
    2218         182 :   cScT.powracpi = powersr(sqrtr(mppi(prec2)), r1);
    2219             : 
    2220         182 :   cScT.cS = cgetg(n0+1, t_VEC);
    2221         182 :   cScT.cT = cgetg(n0+1, t_VEC);
    2222         182 :   for (j=1; j<=n0; j++) gel(cScT.cS,j) = gel(cScT.cT,j) = NULL;
    2223             : 
    2224         182 :   cScT.i0 = i0;
    2225             : 
    2226         182 :   av1 = avma;
    2227         476 :   for (jc = 1; jc <= ncond; jc++)
    2228             :   {
    2229         294 :     const GEN LChar = gel(vChar,jc);
    2230         294 :     const long nChar = lg(LChar)-1, NN = N0[jc];
    2231             : 
    2232         294 :     if (DEBUGLEVEL>1)
    2233           0 :       err_printf("* conductor no %ld/%ld (N = %ld)\n\tInit: ", jc,ncond,NN);
    2234             : 
    2235         294 :     cScT.c1 = gel(C,jc);
    2236         294 :     init_cScT(&cScT, gel(dataCR, LChar[1]), NN, prec2);
    2237         294 :     av2 = avma;
    2238         847 :     for (k = 1; k <= nChar; k++)
    2239             :     {
    2240         553 :       const long t = LChar[k];
    2241         553 :       if (DEBUGLEVEL>1)
    2242           0 :         err_printf("\tcharacter no: %ld (%ld/%ld)\n", t,k,nChar);
    2243             : 
    2244         553 :       if (!isintzero( ch_comp(gel(dataCR, t)) ))
    2245             :       {
    2246         546 :         const long d = degs[t];
    2247         546 :         const GEN dtcr = gel(dataCR, t), z = gel(ch_CHI(dtcr), 2);
    2248         546 :         GEN p1 = gen_0, p2 = gen_0;
    2249         546 :         long c = 0;
    2250         546 :         int **matan = ComputeCoeff(gel(dataCR,t), &LIST, NN, d);
    2251     1794135 :         for (n = 1; n <= NN; n++)
    2252     1793589 :           if ((an = EvalCoeff(z, matan[n], d)))
    2253             :           {
    2254      503195 :            get_cS_cT(&cScT, n);
    2255      503195 :            p1 = gadd(p1, gmul(an, gel(cScT.cS,n)));
    2256      503195 :            p2 = gadd(p2, gmul(an, gel(cScT.cT,n)));
    2257      503195 :            if (++c == 256) { gerepileall(av2,2, &p1,&p2); c = 0; }
    2258             :           }
    2259         546 :         gaffect(p1,        gel(S,t));
    2260         546 :         gaffect(gconj(p2), gel(T,t));
    2261         546 :         FreeMat(matan, NN); avma = av2;
    2262             :       }
    2263           7 :       else if (DEBUGLEVEL>1)
    2264           0 :         err_printf("\t  no need to compute this character\n");
    2265             :     }
    2266         294 :     if (DEBUGLEVEL>1) err_printf("\n");
    2267         294 :     avma = av1;
    2268             :   }
    2269         182 :   clear_cScT(&cScT, n0);
    2270         182 :   avma = av;
    2271         182 : }
    2272             : 
    2273             : static void
    2274         399 : GetST(GEN bnr, GEN *pS, GEN *pT, GEN dataCR, GEN vChar, long prec)
    2275             : {
    2276         399 :   const long cl = lg(dataCR) - 1;
    2277         399 :   GEN S, T, nf  = checknf(bnr);
    2278             :   long j;
    2279             : 
    2280             :   /* allocate memory for answer */
    2281         399 :   *pS = S = cgetg(cl+1, t_VEC);
    2282         399 :   *pT = T = cgetg(cl+1, t_VEC);
    2283        1750 :   for (j = 1; j <= cl; j++)
    2284             :   {
    2285        1351 :     gel(S,j) = cgetc(prec);
    2286        1351 :     gel(T,j) = cgetc(prec);
    2287             :   }
    2288         399 :   if (nf_get_degree(nf) == 2)
    2289         252 :     QuadGetST(bnr, pS, pT, dataCR, vChar, prec);
    2290             :   else
    2291         147 :     GetST0(bnr, pS, pT, dataCR, vChar, prec);
    2292         399 : }
    2293             : 
    2294             : /*******************************************************************/
    2295             : /*                                                                 */
    2296             : /*     Class fields of real quadratic fields using Stark units     */
    2297             : /*                                                                 */
    2298             : /*******************************************************************/
    2299             : /* compute the Hilbert class field using genus class field theory when
    2300             :    the exponent of the class group is exactly 2 (trivial group not covered) */
    2301             : /* Cf Herz, Construction of class fields, LNM 21, Theorem 1 (VII-6) */
    2302             : static GEN
    2303          14 : GenusFieldQuadReal(GEN disc)
    2304             : {
    2305          14 :   long i, i0 = 0, l;
    2306          14 :   pari_sp av = avma;
    2307          14 :   GEN T = NULL, p0 = NULL, P;
    2308             : 
    2309          14 :   P = gel(Z_factor(disc), 1);
    2310          14 :   l = lg(P);
    2311          42 :   for (i = 1; i < l; i++)
    2312             :   {
    2313          35 :     GEN p = gel(P,i);
    2314          35 :     if (mod4(p) == 3) { p0 = p; i0 = i; break; }
    2315             :   }
    2316          14 :   l--; /* remove last prime */
    2317          14 :   if (i0 == l) l--; /* ... remove p0 and last prime */
    2318          49 :   for (i = 1; i < l; i++)
    2319             :   {
    2320          35 :     GEN p = gel(P,i), d, t;
    2321          35 :     if (i == i0) continue;
    2322          28 :     if (absequaliu(p, 2))
    2323          14 :       switch (mod32(disc))
    2324             :       {
    2325          14 :         case  8: d = gen_2; break;
    2326           0 :         case 24: d = shifti(p0, 1); break;
    2327           0 :         default: d = p0; break;
    2328             :       }
    2329             :     else
    2330          14 :       d = (mod4(p) == 1)? p: mulii(p0, p);
    2331          28 :     t = mkpoln(3, gen_1, gen_0, negi(d)); /* x^2 - d */
    2332          28 :     T = T? ZX_compositum_disjoint(T, t): t;
    2333             :   }
    2334          14 :   return gerepileupto(av, polredbest(T, 0));
    2335             : }
    2336             : static GEN
    2337         406 : GenusFieldQuadImag(GEN disc)
    2338             : {
    2339             :   long i, l;
    2340         406 :   pari_sp av = avma;
    2341         406 :   GEN T = NULL, P;
    2342             : 
    2343         406 :   P = gel(absZ_factor(disc), 1);
    2344         406 :   l = lg(P);
    2345         406 :   l--; /* remove last prime */
    2346        1183 :   for (i = 1; i < l; i++)
    2347             :   {
    2348         777 :     GEN p = gel(P,i), d, t;
    2349         777 :     if (absequaliu(p, 2))
    2350         231 :       switch (mod32(disc))
    2351             :       {
    2352          56 :         case 24: d = gen_2; break;  /* disc = 8 mod 32 */
    2353          42 :         case  8: d = gen_m2; break; /* disc =-8 mod 32 */
    2354         133 :         default: d = gen_m1; break;
    2355             :       }
    2356             :     else
    2357         546 :       d = (mod4(p) == 1)? p: negi(p);
    2358         777 :     t = mkpoln(3, gen_1, gen_0, negi(d)); /* x^2 - d */
    2359         777 :     T = T? ZX_compositum_disjoint(T, t): t;
    2360             :   }
    2361         406 :   return gerepileupto(av, polredbest(T, 0));
    2362             : }
    2363             : 
    2364             : /* if flag != 0, computes a fast and crude approximation of the result */
    2365             : static GEN
    2366         644 : AllStark(GEN data,  GEN nf,  long flag,  long newprec)
    2367             : {
    2368         644 :   const long BND = 300;
    2369         644 :   long cl, i, j, cpt = 0, N, h, v, n, r1, r2, den;
    2370             :   pari_sp av, av2;
    2371             :   int **matan;
    2372         644 :   GEN bnr = gel(data,1), p1, p2, S, T, polrelnum, polrel, Lp, W, veczeta;
    2373             :   GEN vChar, degs, C, dataCR, cond1, L1, an;
    2374             :   LISTray LIST;
    2375             :   pari_timer ti;
    2376             : 
    2377         644 :   nf_get_sign(nf, &r1,&r2);
    2378         644 :   N     = nf_get_degree(nf);
    2379         644 :   cond1 = gel(bnr_get_mod(bnr), 2);
    2380         644 :   dataCR = gel(data,5);
    2381         644 :   vChar = sortChars(dataCR);
    2382             : 
    2383         644 :   v = 1;
    2384         644 :   while (gequal1(gel(cond1,v))) v++;
    2385             : 
    2386         644 :   cl = lg(dataCR)-1;
    2387         644 :   degs = GetDeg(dataCR);
    2388         644 :   h  = itos(ZM_det_triangular(gel(data,2))) >> 1;
    2389             : 
    2390             : LABDOUB:
    2391         665 :   if (DEBUGLEVEL) timer_start(&ti);
    2392         665 :   av = avma;
    2393             : 
    2394             :   /* characters with rank > 1 should not be computed */
    2395        2842 :   for (i = 1; i <= cl; i++)
    2396             :   {
    2397        2177 :     GEN chi = gel(dataCR, i);
    2398        2177 :     if (L_vanishes_at_0(chi)) ch_comp(chi) = gen_0;
    2399             :   }
    2400             : 
    2401         665 :   W = ComputeAllArtinNumbers(dataCR, vChar, (flag >= 0), newprec);
    2402         665 :   if (DEBUGLEVEL) timer_printf(&ti,"Compute W");
    2403         665 :   Lp = cgetg(cl + 1, t_VEC);
    2404         665 :   if (!flag)
    2405             :   {
    2406         343 :     GetST(bnr, &S, &T, dataCR, vChar, newprec);
    2407         343 :     if (DEBUGLEVEL) timer_printf(&ti, "S&T");
    2408        1470 :     for (i = 1; i <= cl; i++)
    2409             :     {
    2410        1127 :       GEN chi = gel(dataCR, i), v = gen_0;
    2411        1127 :       if (!isintzero( ch_comp(chi) ))
    2412        1113 :         v = gel(GetValue(chi, gel(W,i), gel(S,i), gel(T,i), 2, newprec), 2);
    2413        1127 :       gel(Lp, i) = v;
    2414             :     }
    2415             :   }
    2416             :   else
    2417             :   { /* compute a crude approximation of the result */
    2418         322 :     C = cgetg(cl + 1, t_VEC);
    2419         322 :     for (i = 1; i <= cl; i++) gel(C,i) = ch_C(gel(dataCR, i));
    2420         322 :     n = zeta_get_N0(vecmax(C), zeta_get_limx(r1, r2, prec2nbits(newprec)));
    2421         322 :     if (n > BND) n = BND;
    2422         322 :     if (DEBUGLEVEL) err_printf("N0 in QuickPol: %ld \n", n);
    2423         322 :     InitPrimes(bnr, n, &LIST);
    2424             : 
    2425         322 :     L1 = cgetg(cl+1, t_VEC);
    2426             :     /* use L(1) = sum (an / n) */
    2427        1372 :     for (i = 1; i <= cl; i++)
    2428             :     {
    2429        1050 :       GEN dtcr = gel(dataCR,i);
    2430        1050 :       matan = ComputeCoeff(dtcr, &LIST, n, degs[i]);
    2431        1050 :       av2 = avma;
    2432        1050 :       p1 = real_0(newprec); p2 = gel(ch_CHI(dtcr), 2);
    2433      310366 :       for (j = 1; j <= n; j++)
    2434      309316 :         if ( (an = EvalCoeff(p2, matan[j], degs[i])) )
    2435      108920 :           p1 = gadd(p1, gdivgs(an, j));
    2436        1050 :       gel(L1,i) = gerepileupto(av2, p1);
    2437        1050 :       FreeMat(matan, n);
    2438             :     }
    2439         322 :     p1 = gmul2n(powruhalf(mppi(newprec), N-2), 1);
    2440             : 
    2441        1372 :     for (i = 1; i <= cl; i++)
    2442             :     {
    2443             :       long r;
    2444        1050 :       GEN WW, A = ComputeAChi(gel(dataCR,i), &r, 0, newprec);
    2445        1050 :       WW = gmul(gel(C,i), gmul(A, gel(W,i)));
    2446        1050 :       gel(Lp,i) = gdiv(gmul(WW, gconj(gel(L1,i))), p1);
    2447             :     }
    2448             :   }
    2449             : 
    2450         665 :   p1 = ComputeLift(gel(data,4));
    2451             : 
    2452         665 :   den = flag ? h: 2*h;
    2453         665 :   veczeta = cgetg(h + 1, t_VEC);
    2454        4214 :   for (i = 1; i <= h; i++)
    2455             :   {
    2456        3549 :     GEN z = gen_0, sig = gel(p1,i);
    2457       18165 :     for (j = 1; j <= cl; j++)
    2458             :     {
    2459       14616 :       GEN dtcr = gel(dataCR,j), CHI = ch_CHI(dtcr);
    2460       14616 :       GEN t = mulreal(gel(Lp,j), CharEval(CHI, sig));
    2461       14616 :       if (chi_get_deg(CHI) != 2) t = gmul2n(t, 1); /* character not real */
    2462       14616 :       z = gadd(z, t);
    2463             :     }
    2464        3549 :     gel(veczeta,i) = gdivgs(z, den);
    2465             :   }
    2466        4214 :   for (j = 1; j <= h; j++)
    2467        3549 :     gel(veczeta,j) = gmul2n(gcosh(gel(veczeta,j), newprec), 1);
    2468         665 :   polrelnum = roots_to_pol(veczeta, 0);
    2469         665 :   if (DEBUGLEVEL)
    2470             :   {
    2471           0 :     if (DEBUGLEVEL>1) {
    2472           0 :       err_printf("polrelnum = %Ps\n", polrelnum);
    2473           0 :       err_printf("zetavalues = %Ps\n", veczeta);
    2474           0 :       if (!flag)
    2475           0 :         err_printf("Checking the square-root of the Stark unit...\n");
    2476             :     }
    2477           0 :     timer_printf(&ti, "Compute %s", flag? "quickpol": "polrelnum");
    2478             :   }
    2479             : 
    2480         665 :   if (flag)
    2481         322 :     return gerepilecopy(av, polrelnum);
    2482             : 
    2483             :   /* try to recognize this polynomial */
    2484         343 :   polrel = RecCoeff(nf, polrelnum, v, newprec);
    2485         343 :   if (!polrel)
    2486             :   {
    2487        1666 :     for (j = 1; j <= h; j++)
    2488        1428 :       gel(veczeta,j) = gsubgs(gsqr(gel(veczeta,j)), 2);
    2489         238 :     polrelnum = roots_to_pol(veczeta, 0);
    2490         238 :     if (DEBUGLEVEL)
    2491             :     {
    2492           0 :       if (DEBUGLEVEL>1) {
    2493           0 :         err_printf("It's not a square...\n");
    2494           0 :         err_printf("polrelnum = %Ps\n", polrelnum);
    2495             :       }
    2496           0 :       timer_printf(&ti, "Compute polrelnum");
    2497             :     }
    2498         238 :     polrel = RecCoeff(nf, polrelnum, v, newprec);
    2499             :   }
    2500         343 :   if (!polrel) /* FAILED */
    2501             :   {
    2502          21 :     const long EXTRA_BITS = 64;
    2503             :     long incr_pr;
    2504          21 :     if (++cpt >= 3) pari_err_PREC( "stark (computation impossible)");
    2505             :     /* estimate needed precision */
    2506          21 :     incr_pr = prec2nbits(gprecision(polrelnum))- gexpo(polrelnum);
    2507          21 :     if (incr_pr < 0) incr_pr = -incr_pr + EXTRA_BITS;
    2508          21 :     newprec += nbits2extraprec(maxss(3*EXTRA_BITS, cpt*incr_pr));
    2509          21 :     if (DEBUGLEVEL) pari_warn(warnprec, "AllStark", newprec);
    2510             : 
    2511          21 :     nf = nfnewprec_shallow(nf, newprec);
    2512          21 :     dataCR = CharNewPrec(dataCR, nf, newprec);
    2513             : 
    2514          21 :     gerepileall(av, 2, &nf, &dataCR);
    2515          21 :     goto LABDOUB;
    2516             :   }
    2517             : 
    2518         322 :   if (DEBUGLEVEL) {
    2519           0 :     if (DEBUGLEVEL>1) err_printf("polrel = %Ps\n", polrel);
    2520           0 :     timer_printf(&ti, "Recpolnum");
    2521             :   }
    2522         322 :   return gerepilecopy(av, polrel);
    2523             : }
    2524             : 
    2525             : /********************************************************************/
    2526             : /*                        Main functions                            */
    2527             : /********************************************************************/
    2528             : 
    2529             : static GEN
    2530         259 : get_subgroup(GEN H, GEN cyc, const char *s)
    2531             : {
    2532         259 :   if (!H || gequal0(H)) return diagonal_shallow(cyc);
    2533          21 :   if (typ(H) != t_MAT) pari_err_TYPE(stack_strcat(s," [subgroup]"), H);
    2534          14 :   RgM_check_ZM(H, s);
    2535          14 :   return ZM_hnfmodid(H, cyc);
    2536             : }
    2537             : 
    2538             : GEN
    2539         203 : bnrstark(GEN bnr, GEN subgrp, long prec)
    2540             : {
    2541             :   long N, newprec;
    2542         203 :   pari_sp av = avma;
    2543             :   GEN bnf, p1, cycbnr, nf, data, dtQ;
    2544             : 
    2545             :   /* check the bnr */
    2546         203 :   checkbnr(bnr);
    2547         203 :   bnf = checkbnf(bnr);
    2548         203 :   nf  = bnf_get_nf(bnf);
    2549         203 :   N   = nf_get_degree(nf);
    2550         203 :   if (N == 1) return galoissubcyclo(bnr, subgrp, 0, 0);
    2551             : 
    2552             :   /* check the bnf */
    2553         203 :   if (!nf_get_varn(nf))
    2554           0 :     pari_err_PRIORITY("bnrstark", nf_get_pol(nf), "=", 0);
    2555         203 :   if (nf_get_r2(nf)) pari_err_DOMAIN("bnrstark", "r2", "!=", gen_0, nf);
    2556         196 :   subgrp = get_subgroup(subgrp,bnr_get_cyc(bnr),"bnrstark");
    2557             : 
    2558             :   /* compute bnr(conductor) */
    2559         196 :   p1     = bnrconductor_i(bnr, subgrp, 2);
    2560         196 :   bnr    = gel(p1,2); cycbnr = bnr_get_cyc(bnr);
    2561         196 :   subgrp = gel(p1,3);
    2562         196 :   if (gequal1( ZM_det_triangular(subgrp) )) { avma = av; return pol_x(0); }
    2563             : 
    2564             :   /* check the class field */
    2565         196 :   if (!gequal0(gel(bnr_get_mod(bnr), 2)))
    2566           7 :     pari_err_DOMAIN("bnrstark", "r2(class field)", "!=", gen_0, bnr);
    2567             : 
    2568             :   /* find a suitable extension N */
    2569         189 :   dtQ = InitQuotient(subgrp);
    2570         189 :   data  = FindModulus(bnr, dtQ, &newprec);
    2571         189 :   if (!data)
    2572             :   {
    2573           0 :     GEN vec, H, cyc = gel(dtQ,2), U = gel(dtQ,3), M = RgM_inv(U);
    2574           0 :     long i, j = 1, l = lg(M);
    2575             : 
    2576             :     /* M = indep. generators of Cl_f/subgp, restrict to cyclic components */
    2577           0 :     vec = cgetg(l, t_VEC);
    2578           0 :     for (i = 1; i < l; i++)
    2579             :     {
    2580           0 :       if (is_pm1(gel(cyc,i))) continue;
    2581           0 :       H = ZM_hnfmodid(vecsplice(M,i), cycbnr);
    2582           0 :       gel(vec,j++) = bnrstark(bnr, H, prec);
    2583             :     }
    2584           0 :     setlg(vec, j); return gerepilecopy(av, vec);
    2585             :   }
    2586             : 
    2587         189 :   if (newprec > prec)
    2588             :   {
    2589          27 :     if (DEBUGLEVEL>1) err_printf("new precision: %ld\n", newprec);
    2590          27 :     nf = nfnewprec_shallow(nf, newprec);
    2591             :   }
    2592         189 :   return gerepileupto(av, AllStark(data, nf, 0, newprec));
    2593             : }
    2594             : 
    2595             : /* For each character of Cl(bnr)/subgp, compute L(1, chi) (or equivalently
    2596             :  * the first non-zero term c(chi) of the expansion at s = 0).
    2597             :  * If flag & 1: compute the value at s = 1 (for non-trivial characters),
    2598             :  * else compute the term c(chi) and return [r(chi), c(chi)] where r(chi) is
    2599             :  *   the order of L(s, chi) at s = 0.
    2600             :  * If flag & 2: compute the value of the L-function L_S(s, chi) where S is the
    2601             :  *   set of places dividing the modulus of bnr (and the infinite places),
    2602             :  * else
    2603             :  *   compute the value of the primitive L-function attached to chi,
    2604             :  * If flag & 4: return also the character */
    2605             : GEN
    2606          63 : bnrL1(GEN bnr, GEN subgp, long flag, long prec)
    2607             : {
    2608             :   GEN cyc, L1, allCR, listCR;
    2609             :   GEN indCR, invCR, Qt;
    2610             :   long cl, i, nc;
    2611          63 :   pari_sp av = avma;
    2612             : 
    2613          63 :   checkbnr(bnr);
    2614          63 :   if (flag < 0 || flag > 8) pari_err_FLAG("bnrL1");
    2615             : 
    2616          63 :   cyc  = bnr_get_cyc(bnr);
    2617          63 :   subgp = get_subgroup(subgp, cyc, "bnrL1");
    2618             : 
    2619          56 :   Qt = InitQuotient(subgp);
    2620          56 :   cl = itou(gel(Qt,1));
    2621             : 
    2622             :   /* compute all characters */
    2623          56 :   allCR = EltsOfGroup(cl, gel(Qt,2));
    2624             : 
    2625             :   /* make a list of all non-trivial characters modulo conjugation */
    2626          56 :   listCR = cgetg(cl, t_VEC);
    2627          56 :   indCR = cgetg(cl, t_VECSMALL);
    2628          56 :   invCR = cgetg(cl, t_VECSMALL); nc = 0;
    2629         406 :   for (i = 1; i < cl; i++)
    2630             :   {
    2631             :     /* lift to a character on Cl(bnr) */
    2632         350 :     GEN lchi = LiftChar(Qt, cyc, gel(allCR,i));
    2633         350 :     GEN clchi = charconj(cyc, lchi);
    2634         350 :     long j, a = 0;
    2635        1351 :     for (j = 1; j <= nc; j++)
    2636        1127 :       if (ZV_equal(gmael(listCR, j, 1), clchi)) { a = j; break; }
    2637             : 
    2638         350 :     if (!a)
    2639             :     {
    2640         224 :       nc++;
    2641         224 :       gel(listCR,nc) = mkvec2(lchi, bnrconductorofchar(bnr, lchi));
    2642         224 :       indCR[i]  = nc;
    2643         224 :       invCR[nc] = i;
    2644             :     }
    2645             :     else
    2646         126 :       indCR[i] = -invCR[a];
    2647             : 
    2648         350 :     gel(allCR,i) = lchi;
    2649             :   }
    2650          56 :   settyp(allCR[cl], t_VEC); /* set correct type for trivial character */
    2651             : 
    2652          56 :   setlg(listCR, nc + 1);
    2653          56 :   L1 = cgetg((flag&1)? cl: cl+1, t_VEC);
    2654          56 :   if (nc)
    2655             :   {
    2656          56 :     GEN dataCR = InitChar(bnr, listCR, prec);
    2657          56 :     GEN W, S, T, vChar = sortChars(dataCR);
    2658          56 :     GetST(bnr, &S, &T, dataCR, vChar, prec);
    2659          56 :     W = ComputeAllArtinNumbers(dataCR, vChar, 1, prec);
    2660         406 :     for (i = 1; i < cl; i++)
    2661             :     {
    2662         350 :       long a = indCR[i];
    2663         350 :       if (a > 0)
    2664         224 :         gel(L1,i) = GetValue(gel(dataCR,a), gel(W,a), gel(S,a), gel(T,a),
    2665             :                              flag, prec);
    2666             :       else
    2667         126 :         gel(L1,i) = gconj(gel(L1,-a));
    2668             :     }
    2669             :   }
    2670          56 :   if (!(flag & 1))
    2671          21 :     gel(L1,cl) = GetValue1(bnr, flag & 2, prec);
    2672             :   else
    2673          35 :     cl--;
    2674             : 
    2675          56 :   if (flag & 4) {
    2676          28 :     for (i = 1; i <= cl; i++) gel(L1,i) = mkvec2(gel(allCR,i), gel(L1,i));
    2677             :   }
    2678          56 :   return gerepilecopy(av, L1);
    2679             : }
    2680             : 
    2681             : /*******************************************************************/
    2682             : /*                                                                 */
    2683             : /*       Hilbert and Ray Class field using Stark                   */
    2684             : /*                                                                 */
    2685             : /*******************************************************************/
    2686             : /* P in A[x,y], deg_y P < 2, return P0 and P1 in A[x] such that P = P0 + P1 y */
    2687             : static void
    2688         133 : split_pol_quad(GEN P, GEN *gP0, GEN *gP1)
    2689             : {
    2690         133 :   long i, l = lg(P);
    2691         133 :   GEN P0 = cgetg(l, t_POL), P1 = cgetg(l, t_POL);
    2692         133 :   P0[1] = P1[1] = P[1];
    2693        1211 :   for (i = 2; i < l; i++)
    2694             :   {
    2695        1078 :     GEN c = gel(P,i), c0 = c, c1 = gen_0;
    2696        1078 :     if (typ(c) == t_POL) /* write c = c1 y + c0 */
    2697         945 :       switch(degpol(c))
    2698             :       {
    2699           0 :         case -1: c0 = gen_0; break;
    2700         945 :         default: c1 = gel(c,3); /* fall through */
    2701         945 :         case  0: c0 = gel(c,2); break;
    2702             :       }
    2703        1078 :     gel(P0,i) = c0; gel(P1,i) = c1;
    2704             :   }
    2705         133 :   *gP0 = normalizepol_lg(P0, l);
    2706         133 :   *gP1 = normalizepol_lg(P1, l);
    2707         133 : }
    2708             : 
    2709             : /* k = nf quadratic field, P relative equation of H_k (Hilbert class field)
    2710             :  * return T in Z[X], such that H_k / Q is the compositum of Q[X]/(T) and k */
    2711             : static GEN
    2712         133 : makescind(GEN nf, GEN P)
    2713             : {
    2714         133 :   GEN Pp, p, pol, G, L, a, roo, P0,P1, Ny,Try, nfpol = nf_get_pol(nf);
    2715             :   long i, is_P;
    2716             : 
    2717         133 :   P = lift_shallow(P);
    2718         133 :   split_pol_quad(P, &P0, &P1);
    2719             :   /* P = P0 + y P1, Norm_{k/Q}(P) = P0^2 + Tr y P0P1 + Ny P1^2, irreducible/Q */
    2720         133 :   Ny = gel(nfpol, 2);
    2721         133 :   Try = negi(gel(nfpol, 3));
    2722         133 :   pol = RgX_add(RgX_sqr(P0), RgX_Rg_mul(RgX_sqr(P1), Ny));
    2723         133 :   if (signe(Try)) pol = RgX_add(pol, RgX_Rg_mul(RgX_mul(P0,P1), Try));
    2724             :   /* pol = rnfequation(nf, P); */
    2725         133 :   G = galoisinit(pol, NULL);
    2726         133 :   L = gal_get_group(G);
    2727         133 :   p = gal_get_p(G);
    2728         133 :   a = FpX_quad_root(nfpol, p, 0);
    2729             :   /* P mod a prime \wp above p (which splits) */
    2730         133 :   Pp = FpXY_evalx(P, a, p);
    2731         133 :   roo = gal_get_roots(G);
    2732         133 :   is_P = gequal0( FpX_eval(Pp, remii(gel(roo,1),p), p) );
    2733             :   /* each roo[i] mod p is a root of P or (exclusive) tau(P) mod \wp */
    2734             :   /* record whether roo[1] is a root of P or tau(P) */
    2735             : 
    2736        1078 :   for (i = 1; i < lg(L); i++)
    2737             :   {
    2738        1078 :     GEN perm = gel(L,i);
    2739        1078 :     long k = perm[1]; if (k == 1) continue;
    2740         945 :     k = gequal0( FpX_eval(Pp, remii(gel(roo,k),p), p) );
    2741             :     /* roo[k] is a root of the other polynomial */
    2742         945 :     if (k != is_P)
    2743             :     {
    2744         133 :       long o = perm_order(perm);
    2745         133 :       if (o != 2) perm = perm_pow(perm, o >> 1);
    2746             :       /* perm has order two and doesn't belong to Gal(H_k/k) */
    2747         133 :       return galoisfixedfield(G, perm, 1, varn(P));
    2748             :     }
    2749             :   }
    2750           0 :   pari_err_BUG("makescind");
    2751           0 :   return NULL; /*not reached*/
    2752             : }
    2753             : 
    2754             : /* pbnf = NULL if no bnf is needed, f = NULL may be passed for a trivial
    2755             :  * conductor */
    2756             : static void
    2757         847 : quadray_init(GEN *pD, GEN f, GEN *pbnf, long prec)
    2758             : {
    2759         847 :   GEN D = *pD, nf, bnf = NULL;
    2760         847 :   if (typ(D) == t_INT)
    2761             :   {
    2762             :     int isfund;
    2763         812 :     if (pbnf) {
    2764         252 :       long v = f? gvar(f): NO_VARIABLE;
    2765         252 :       if (v == NO_VARIABLE) v = 1;
    2766         252 :       bnf = Buchall(quadpoly0(D, v), nf_FORCE, prec);
    2767         252 :       nf = bnf_get_nf(bnf);
    2768         252 :       isfund = equalii(D, nf_get_disc(nf));
    2769             :     }
    2770             :     else
    2771         560 :       isfund = Z_isfundamental(D);
    2772         812 :     if (!isfund) pari_err_DOMAIN("quadray", "isfundamental(D)", "=",gen_0, D);
    2773             :   }
    2774             :   else
    2775             :   {
    2776          35 :     bnf = checkbnf(D);
    2777          35 :     nf = bnf_get_nf(bnf);
    2778          35 :     if (nf_get_degree(nf) != 2)
    2779           7 :       pari_err_DOMAIN("quadray", "degree", "!=", gen_2, nf_get_pol(nf));
    2780          28 :     D = nf_get_disc(nf);
    2781             :   }
    2782         833 :   if (pbnf) *pbnf = bnf;
    2783         833 :   *pD = D;
    2784         833 : }
    2785             : 
    2786             : /* compute the polynomial over Q of the Hilbert class field of
    2787             :    Q(sqrt(D)) where D is a positive fundamental discriminant */
    2788             : static GEN
    2789         147 : quadhilbertreal(GEN D, long prec)
    2790             : {
    2791         147 :   pari_sp av = avma;
    2792             :   long newprec;
    2793             :   GEN bnf;
    2794             :   VOLATILE GEN bnr, dtQ, data, nf, cyc, M;
    2795             :   pari_timer ti;
    2796         147 :   if (DEBUGLEVEL) timer_start(&ti);
    2797             : 
    2798             :   (void)&prec; /* prevent longjmp clobbering it */
    2799             :   (void)&bnf;  /* prevent longjmp clobbering it, avoid warning due to
    2800             :                 * quadray_init call : discards qualifiers from pointer type */
    2801         147 :   quadray_init(&D, NULL, &bnf, prec);
    2802         147 :   cyc = bnf_get_cyc(bnf);
    2803         147 :   if (lg(cyc) == 1) { avma = av; return pol_x(0); }
    2804             :   /* if the exponent of the class group is 2, use Genus Theory */
    2805         147 :   if (absequaliu(gel(cyc,1), 2)) return gerepileupto(av, GenusFieldQuadReal(D));
    2806             : 
    2807         133 :   bnr  = Buchray(bnf, gen_1, nf_INIT);
    2808         133 :   M = diagonal_shallow(bnr_get_cyc(bnr));
    2809         133 :   dtQ = InitQuotient(M);
    2810         133 :   nf  = bnf_get_nf(bnf);
    2811             : 
    2812             :   for(;;) {
    2813         133 :     VOLATILE GEN pol = NULL;
    2814         133 :     pari_CATCH(e_PREC) {
    2815           0 :       prec += EXTRA_PREC;
    2816           0 :       if (DEBUGLEVEL) pari_warn(warnprec, "quadhilbertreal", prec);
    2817           0 :       bnr = bnrnewprec_shallow(bnr, prec);
    2818           0 :       bnf = bnr_get_bnf(bnr);
    2819           0 :       nf  = bnf_get_nf(bnf);
    2820             :     } pari_TRY {
    2821             :       /* find the modulus defining N */
    2822             :       pari_timer T;
    2823         133 :       if (DEBUGLEVEL) timer_start(&T);
    2824         133 :       data = FindModulus(bnr, dtQ, &newprec);
    2825         133 :       if (DEBUGLEVEL) timer_printf(&T,"FindModulus");
    2826         133 :       if (!data)
    2827             :       {
    2828           0 :         long i, l = lg(M);
    2829           0 :         GEN vec = cgetg(l, t_VEC);
    2830           0 :         for (i = 1; i < l; i++)
    2831             :         {
    2832           0 :           GEN t = gcoeff(M,i,i);
    2833           0 :           gcoeff(M,i,i) = gen_1;
    2834           0 :           gel(vec,i) = bnrstark(bnr, M, prec);
    2835           0 :           gcoeff(M,i,i) = t;
    2836             :         }
    2837           0 :         return gerepileupto(av, vec);
    2838             :       }
    2839             : 
    2840         133 :       if (newprec > prec)
    2841             :       {
    2842          19 :         if (DEBUGLEVEL>1) err_printf("new precision: %ld\n", newprec);
    2843          19 :         nf = nfnewprec_shallow(nf, newprec);
    2844             :       }
    2845         133 :       pol = AllStark(data, nf, 0, newprec);
    2846         133 :     } pari_ENDCATCH;
    2847         133 :     if (pol) {
    2848         133 :       pol = makescind(nf, pol);
    2849         133 :       return gerepileupto(av, polredbest(pol, 0));
    2850             :     }
    2851           0 :   }
    2852             : }
    2853             : 
    2854             : /*******************************************************************/
    2855             : /*                                                                 */
    2856             : /*       Hilbert and Ray Class field using CM (Schertz)            */
    2857             : /*                                                                 */
    2858             : /*******************************************************************/
    2859             : /* form^2 = 1 ? */
    2860             : static int
    2861         813 : hasexp2(GEN form)
    2862             : {
    2863         813 :   GEN a = gel(form,1), b = gel(form,2), c = gel(form,3);
    2864         813 :   return !signe(b) || absequalii(a,b) || equalii(a,c);
    2865             : }
    2866             : static int
    2867        1323 : uhasexp2(GEN form)
    2868             : {
    2869        1323 :   long a = form[1], b = form[2], c = form[3];
    2870        1323 :   return !b || a == labs(b) || a == c;
    2871             : }
    2872             : 
    2873             : GEN
    2874         455 : qfbforms(GEN D)
    2875             : {
    2876         455 :   ulong d = itou(D), dover3 = d/3, t, b2, a, b, c, h;
    2877         455 :   GEN L = cgetg((long)(sqrt((double)d) * log2(d)), t_VEC);
    2878         455 :   b2 = b = (d&1); h = 0;
    2879         455 :   if (!b) /* b = 0 treated separately to avoid special cases */
    2880             :   {
    2881         252 :     t = d >> 2; /* (b^2 - D) / 4*/
    2882        2954 :     for (a=1; a*a<=t; a++)
    2883        2702 :       if (c = t/a, t == c*a) gel(L,++h) = mkvecsmall3(a,0,c);
    2884         252 :     b = 2; b2 = 4;
    2885             :   }
    2886             :   /* now b > 0, b = D (mod 2) */
    2887        8078 :   for ( ; b2 <= dover3; b += 2, b2 = b*b)
    2888             :   {
    2889        7623 :     t = (b2 + d) >> 2; /* (b^2 - D) / 4*/
    2890             :     /* b = a */
    2891        7623 :     if (c = t/b, t == c*b) gel(L,++h) = mkvecsmall3(b,b,c);
    2892             :     /* b < a < c */
    2893     1912029 :     for (a = b+1; a*a < t; a++)
    2894     1904406 :       if (c = t/a, t == c*a)
    2895             :       {
    2896        1057 :         gel(L,++h) = mkvecsmall3(a, b,c);
    2897        1057 :         gel(L,++h) = mkvecsmall3(a,-b,c);
    2898             :       }
    2899             :     /* a = c */
    2900        7623 :     if (a * a == t) gel(L,++h) = mkvecsmall3(a,b,a);
    2901             :   }
    2902         455 :   setlg(L,h+1); return L;
    2903             : }
    2904             : 
    2905             : /* gcd(n, 24) */
    2906             : static long
    2907         813 : GCD24(long n)
    2908             : {
    2909         813 :   switch(n % 24)
    2910             :   {
    2911          35 :     case 0: return 24;
    2912          35 :     case 1: return 1;
    2913          28 :     case 2: return 2;
    2914           0 :     case 3: return 3;
    2915         119 :     case 4: return 4;
    2916           0 :     case 5: return 1;
    2917         105 :     case 6: return 6;
    2918           0 :     case 7: return 1;
    2919           0 :     case 8: return 8;
    2920           0 :     case 9: return 3;
    2921          91 :     case 10: return 2;
    2922           0 :     case 11: return 1;
    2923         119 :     case 12: return 12;
    2924           0 :     case 13: return 1;
    2925           0 :     case 14: return 2;
    2926           0 :     case 15: return 3;
    2927          91 :     case 16: return 8;
    2928           0 :     case 17: return 1;
    2929          92 :     case 18: return 6;
    2930           0 :     case 19: return 1;
    2931           0 :     case 20: return 4;
    2932           0 :     case 21: return 3;
    2933          98 :     case 22: return 2;
    2934           0 :     case 23: return 1;
    2935           0 :     default: return 0;
    2936             :   }
    2937             : }
    2938             : 
    2939             : struct gpq_data {
    2940             :   long p, q;
    2941             :   GEN sqd; /* sqrt(D), t_REAL */
    2942             :   GEN u, D;
    2943             :   GEN pq, pq2; /* p*q, 2*p*q */
    2944             :   GEN qfpq ; /* class of \P * \Q */
    2945             : };
    2946             : 
    2947             : /* find P and Q two non principal prime ideals (above p <= q) such that
    2948             :  *   cl(P) = cl(Q) if P,Q have order 2 in Cl(K).
    2949             :  *   Ensure that e = 24 / gcd(24, (p-1)(q-1)) = 1 */
    2950             : /* D t_INT, discriminant */
    2951             : static void
    2952          49 : init_pq(GEN D, struct gpq_data *T)
    2953             : {
    2954          49 :   const long Np = 6547; /* N.B. primepi(50000) = 5133 */
    2955          49 :   const ulong maxq = 50000;
    2956          49 :   GEN listp = cgetg(Np + 1, t_VECSMALL); /* primes p */
    2957          49 :   GEN listP = cgetg(Np + 1, t_VEC); /* primeform(p) if of order 2, else NULL */
    2958          49 :   GEN gcd24 = cgetg(Np + 1, t_VECSMALL); /* gcd(p-1, 24) */
    2959             :   forprime_t S;
    2960          49 :   long l = 1;
    2961          49 :   double best = 0.;
    2962             :   ulong q;
    2963             : 
    2964          49 :   u_forprime_init(&S, 2, ULONG_MAX);
    2965          49 :   T->D = D;
    2966          49 :   T->p = T->q = 0;
    2967             :   for(;;)
    2968             :   {
    2969             :     GEN Q;
    2970             :     long i, gcdq, mod;
    2971             :     int order2, store;
    2972             :     double t;
    2973             : 
    2974        1826 :     q = u_forprime_next(&S);
    2975        1826 :     if (best > 0 && q >= maxq)
    2976             :     {
    2977           0 :       if (DEBUGLEVEL)
    2978           0 :         pari_warn(warner,"possibly suboptimal (p,q) for D = %Ps", D);
    2979           0 :       break;
    2980             :     }
    2981        1826 :     if (kroiu(D, q) < 0) continue; /* inert */
    2982         890 :     Q = redimag(primeform_u(D, q));
    2983         890 :     if (is_pm1(gel(Q,1))) continue; /* Q | q is principal */
    2984             : 
    2985         813 :     store = 1;
    2986         813 :     order2 = hasexp2(Q);
    2987         813 :     gcd24[l] = gcdq = GCD24(q-1);
    2988         813 :     mod = 24 / gcdq; /* mod must divide p-1 otherwise e > 1 */
    2989         813 :     listp[l] = q;
    2990         813 :     gel(listP,l) = order2 ? Q : NULL;
    2991         813 :     t = (q+1)/(double)(q-1);
    2992        2129 :     for (i = 1; i < l; i++) /* try all (p, q), p < q in listp */
    2993             :     {
    2994        1660 :       long p = listp[i], gcdp = gcd24[i];
    2995             :       double b;
    2996             :       /* P,Q order 2 => cl(Q) = cl(P) */
    2997        1660 :       if (order2 && gel(listP,i) && !gequal(gel(listP,i), Q)) continue;
    2998        1653 :       if (gcdp % gcdq == 0) store = 0; /* already a better one in the list */
    2999        1653 :       if ((p-1) % mod) continue;
    3000             : 
    3001         344 :       b = (t*(p+1)) / (p-1); /* (p+1)(q+1) / (p-1)(q-1) */
    3002         344 :       if (b > best) {
    3003          98 :         store = 0; /* (p,q) always better than (q,r) for r >= q */
    3004          98 :         best = b; T->q = q; T->p = p;
    3005          98 :         if (DEBUGLEVEL>2) err_printf("p,q = %ld,%ld\n", p, q);
    3006             :       }
    3007             :       /* won't improve with this q as largest member */
    3008         344 :       if (best > 0) break;
    3009             :     }
    3010             :     /* if !store or (q,r) won't improve on current best pair, forget that q */
    3011         813 :     if (store && t*t > best)
    3012         119 :       if (++l >= Np) pari_err_BUG("quadhilbert (not enough primes)");
    3013         813 :     if (!best) /* (p,q) with p < q always better than (q,q) */
    3014             :     { /* try (q,q) */
    3015         140 :       if (gcdq >= 12 && umodiu(D, q)) /* e = 1 and unramified */
    3016             :       {
    3017           7 :         double b = (t*q) / (q-1); /* q(q+1) / (q-1)^2 */
    3018           7 :         if (b > best) {
    3019           7 :           best = b; T->q = T->p = q;
    3020           7 :           if (DEBUGLEVEL>2) err_printf("p,q = %ld,%ld\n", q, q);
    3021             :         }
    3022             :       }
    3023             :     }
    3024             :     /* If (p1+1)(q+1) / (p1-1)(q-1) <= best, we can no longer improve
    3025             :      * even with best p : stop */
    3026         813 :     if ((listp[1]+1)*t <= (listp[1]-1)*best) break;
    3027        1777 :   }
    3028          49 :   if (DEBUGLEVEL>1)
    3029           0 :     err_printf("(p, q) = %ld, %ld; gain = %f\n", T->p, T->q, 12*best);
    3030          49 : }
    3031             : 
    3032             : static GEN
    3033        4102 : gpq(GEN form, struct gpq_data *T)
    3034             : {
    3035        4102 :   pari_sp av = avma;
    3036        4102 :   long a = form[1], b = form[2], c = form[3];
    3037        4102 :   long p = T->p, q = T->q;
    3038             :   GEN form2, w, z;
    3039        4102 :   int fl, real = 0;
    3040             : 
    3041        4102 :   form2 = qficomp(T->qfpq, mkvec3s(a, -b, c));
    3042             :   /* form2 and form yield complex conjugate roots : only compute for the
    3043             :    * lexicographically smallest of the 2 */
    3044        4102 :   fl = cmpis(gel(form2,1), a);
    3045        4102 :   if (fl <= 0)
    3046             :   {
    3047        2156 :     if (fl < 0) return NULL;
    3048         210 :     fl = cmpis(gel(form2,2), b);
    3049         210 :     if (fl <= 0)
    3050             :     {
    3051         147 :       if (fl < 0) return NULL;
    3052             :       /* form == form2 : real root */
    3053          84 :       real = 1;
    3054             :     }
    3055             :   }
    3056             : 
    3057        2093 :   if (p == 2) { /* (a,b,c) = (1,1,0) mod 2 ? */
    3058         203 :     if (a % q == 0 && (a & b & 1) && !(c & 1))
    3059             :     { /* apply S : make sure that (a,b,c) represents odd values */
    3060           0 :       lswap(a,c); b = -b;
    3061             :     }
    3062             :   }
    3063        2093 :   if (a % p == 0 || a % q == 0)
    3064             :   { /* apply T^k, look for c' = a k^2 + b k + c coprime to N */
    3065        1092 :     while (c % p == 0 || c % q == 0)
    3066             :     {
    3067          98 :       c += a + b;
    3068          98 :       b += a << 1;
    3069             :     }
    3070         497 :     lswap(a, c); b = -b; /* apply S */
    3071             :   }
    3072             :   /* now (a,b,c) ~ form and (a,pq) = 1 */
    3073             : 
    3074             :   /* gcd(2a, u) = 2,  w = u mod 2pq, -b mod 2a */
    3075        2093 :   w = Z_chinese(T->u, stoi(-b), T->pq2, utoipos(a << 1));
    3076        2093 :   z = double_eta_quotient(utoipos(a), w, T->D, T->p, T->q, T->pq, T->sqd);
    3077        2093 :   if (real && typ(z) == t_COMPLEX) z = gcopy(gel(z, 1));
    3078        2093 :   return gerepileupto(av, z);
    3079             : }
    3080             : 
    3081             : /* returns an equation for the Hilbert class field of Q(sqrt(D)), D < 0
    3082             :  * fundamental discriminant */
    3083             : static GEN
    3084         462 : quadhilbertimag(GEN D)
    3085             : {
    3086             :   GEN L, P, Pi, Pr, qfp, u;
    3087         462 :   pari_sp av = avma;
    3088             :   long h, i, prec;
    3089             :   struct gpq_data T;
    3090             :   pari_timer ti;
    3091             : 
    3092         462 :   if (DEBUGLEVEL>1) timer_start(&ti);
    3093         462 :   if (lgefint(D) == 3)
    3094         462 :     switch (D[2]) { /* = |D|; special cases where e > 1 */
    3095             :       case 3:
    3096             :       case 4:
    3097             :       case 7:
    3098             :       case 8:
    3099             :       case 11:
    3100             :       case 19:
    3101             :       case 43:
    3102             :       case 67:
    3103           7 :       case 163: return pol_x(0);
    3104             :     }
    3105         455 :   L = qfbforms(D);
    3106         455 :   h = lg(L)-1;
    3107         455 :   if ((1L << vals(h)) == h) /* power of 2 */
    3108             :   { /* check whether > |Cl|/2 elements have order <= 2 ==> 2-elementary */
    3109         413 :     long lim = (h>>1) + 1;
    3110        1729 :     for (i=1; i <= lim; i++)
    3111        1323 :       if (!uhasexp2(gel(L,i))) break;
    3112         413 :     if (i > lim) return GenusFieldQuadImag(D);
    3113             :   }
    3114          49 :   if (DEBUGLEVEL>1) timer_printf(&ti,"class number = %ld",h);
    3115          49 :   init_pq(D, &T);
    3116          49 :   qfp = primeform_u(D, T.p);
    3117          49 :   T.pq =  muluu(T.p, T.q);
    3118          49 :   T.pq2 = shifti(T.pq,1);
    3119          49 :   if (T.p == T.q)
    3120             :   {
    3121           0 :     GEN qfbp2 = qficompraw(qfp, qfp);
    3122           0 :     u = gel(qfbp2,2);
    3123           0 :     T.u = modii(u, T.pq2);
    3124           0 :     T.qfpq = redimag(qfbp2);
    3125             :   }
    3126             :   else
    3127             :   {
    3128          49 :     GEN qfq = primeform_u(D, T.q), bp = gel(qfp,2), bq = gel(qfq,2);
    3129          49 :     T.u = Z_chinese(bp, bq, utoipos(T.p << 1), utoipos(T.q << 1));
    3130             :     /* T.u = bp (mod 2p), T.u = bq (mod 2q) */
    3131          49 :     T.qfpq = qficomp(qfp, qfq);
    3132             :   }
    3133             :   /* u modulo 2pq */
    3134          49 :   prec = LOWDEFAULTPREC;
    3135          49 :   Pr = cgetg(h+1,t_VEC);
    3136          49 :   Pi = cgetg(h+1,t_VEC);
    3137             :   for(;;)
    3138             :   {
    3139          63 :     long ex, exmax = 0, r1 = 0, r2 = 0;
    3140          63 :     pari_sp av0 = avma;
    3141          63 :     T.sqd = sqrtr_abs(itor(D, prec));
    3142        4165 :     for (i=1; i<=h; i++)
    3143             :     {
    3144        4102 :       GEN s = gpq(gel(L,i), &T);
    3145        4102 :       if (DEBUGLEVEL>3) err_printf("%ld ", i);
    3146        4102 :       if (!s) continue;
    3147        2093 :       if (typ(s) != t_COMPLEX) gel(Pr, ++r1) = s; /* real root */
    3148        2009 :       else                     gel(Pi, ++r2) = s;
    3149        2093 :       ex = gexpo(s); if (ex > 0) exmax += ex;
    3150             :     }
    3151          63 :     if (DEBUGLEVEL>1) timer_printf(&ti,"roots");
    3152          63 :     setlg(Pr, r1+1);
    3153          63 :     setlg(Pi, r2+1);
    3154          63 :     P = roots_to_pol_r1(shallowconcat(Pr,Pi), 0, r1);
    3155          63 :     P = grndtoi(P,&exmax);
    3156          63 :     if (DEBUGLEVEL>1) timer_printf(&ti,"product, error bits = %ld",exmax);
    3157          63 :     if (exmax <= -10) break;
    3158          14 :     avma = av0; prec += nbits2extraprec(prec2nbits(DEFAULTPREC)+exmax);
    3159          14 :     if (DEBUGLEVEL) pari_warn(warnprec,"quadhilbertimag",prec);
    3160          14 :   }
    3161          49 :   return gerepileupto(av,P);
    3162             : }
    3163             : 
    3164             : GEN
    3165         574 : quadhilbert(GEN D, long prec)
    3166             : {
    3167         574 :   GEN d = D;
    3168         574 :   quadray_init(&d, NULL, NULL, 0);
    3169        1120 :   return (signe(d)>0)? quadhilbertreal(D,prec)
    3170         560 :                      : quadhilbertimag(d);
    3171             : }
    3172             : 
    3173             : /* return a vector of all roots of 1 in bnf [not necessarily quadratic] */
    3174             : static GEN
    3175          70 : getallrootsof1(GEN bnf)
    3176             : {
    3177          70 :   GEN T, u, nf = bnf_get_nf(bnf), tu;
    3178          70 :   long i, n = bnf_get_tuN(bnf);
    3179             : 
    3180          70 :   if (n == 2) {
    3181          56 :     long N = nf_get_degree(nf);
    3182          56 :     return mkvec2(scalarcol_shallow(gen_m1, N),
    3183             :                   scalarcol_shallow(gen_1, N));
    3184             :   }
    3185          14 :   tu = poltobasis(nf, bnf_get_tuU(bnf));
    3186          14 :   T = zk_multable(nf, tu);
    3187          14 :   u = cgetg(n+1, t_VEC); gel(u,1) = tu;
    3188          14 :   for (i=2; i <= n; i++) gel(u,i) = ZM_ZC_mul(T, gel(u,i-1));
    3189          14 :   return u;
    3190             : }
    3191             : /* assume bnr has the right conductor */
    3192             : static GEN
    3193          70 : get_lambda(GEN bnr)
    3194             : {
    3195          70 :   GEN bnf = bnr_get_bnf(bnr), nf = bnf_get_nf(bnf), pol = nf_get_pol(nf);
    3196          70 :   GEN f = gel(bnr_get_mod(bnr), 1), labas, lamodf, u;
    3197          70 :   long a, b, f2, i, lu, v = varn(pol);
    3198             : 
    3199          70 :   f2 = 2 * itos(gcoeff(f,1,1));
    3200          70 :   u = getallrootsof1(bnf); lu = lg(u);
    3201         238 :   for (i=1; i<lu; i++)
    3202         168 :     gel(u,i) = ZC_hnfrem(gel(u,i), f); /* roots of 1, mod f */
    3203          70 :   if (DEBUGLEVEL>1)
    3204           0 :     err_printf("quadray: looking for [a,b] != unit mod 2f\n[a,b] = ");
    3205         168 :   for (a=0; a<f2; a++)
    3206        2576 :     for (b=0; b<f2; b++)
    3207             :     {
    3208        2478 :       GEN la = deg1pol_shallow(stoi(a), stoi(b), v); /* ax + b */
    3209        2478 :       if (umodiu(gnorm(mkpolmod(la, pol)), f2) != 1) continue;
    3210         224 :       if (DEBUGLEVEL>1) err_printf("[%ld,%ld] ",a,b);
    3211             : 
    3212         224 :       labas = poltobasis(nf, la);
    3213         224 :       lamodf = ZC_hnfrem(labas, f);
    3214         469 :       for (i=1; i<lu; i++)
    3215         399 :         if (ZV_equal(lamodf, gel(u,i))) break;
    3216         224 :       if (i < lu) continue; /* la = unit mod f */
    3217          70 :       if (DEBUGLEVEL)
    3218             :       {
    3219           0 :         if (DEBUGLEVEL>1) err_printf("\n");
    3220           0 :         err_printf("lambda = %Ps\n",la);
    3221             :       }
    3222          70 :       return labas;
    3223             :     }
    3224           0 :   pari_err_BUG("get_lambda");
    3225           0 :   return NULL;
    3226             : }
    3227             : 
    3228             : static GEN
    3229        8778 : to_approx(GEN nf, GEN a)
    3230             : {
    3231        8778 :   GEN M = nf_get_M(nf);
    3232        8778 :   return gadd(gel(a,1), gmul(gcoeff(M,1,2),gel(a,2)));
    3233             : }
    3234             : /* Z-basis for a (over C) */
    3235             : static GEN
    3236        4354 : get_om(GEN nf, GEN a) {
    3237        4354 :   return mkvec2(to_approx(nf,gel(a,2)),
    3238        4354 :                 to_approx(nf,gel(a,1)));
    3239             : }
    3240             : 
    3241             : /* Compute all elts in class group G = [|G|,c,g], c=cyclic factors, g=gens.
    3242             :  * Set list[j + 1] = g1^e1...gk^ek where j is the integer
    3243             :  *   ek + ck [ e(k-1) + c(k-1) [... + c2 [e1]]...] */
    3244             : static GEN
    3245          70 : getallelts(GEN bnr)
    3246             : {
    3247             :   GEN nf, C, c, g, list, pows, gk;
    3248             :   long lc, i, j, no;
    3249             : 
    3250          70 :   nf = bnr_get_nf(bnr);
    3251          70 :   no = itos( bnr_get_no(bnr) );
    3252          70 :   c = bnr_get_cyc(bnr);
    3253          70 :   g = bnr_get_gen_nocheck(bnr); lc = lg(c)-1;
    3254          70 :   list = cgetg(no+1,t_VEC);
    3255          70 :   gel(list,1) = matid(nf_get_degree(nf)); /* (1) */
    3256          70 :   if (!no) return list;
    3257             : 
    3258          70 :   pows = cgetg(lc+1,t_VEC);
    3259          70 :   c = leafcopy(c); settyp(c, t_VECSMALL);
    3260         140 :   for (i=1; i<=lc; i++)
    3261             :   {
    3262          70 :     long k = itos(gel(c,i));
    3263          70 :     c[i] = k;
    3264          70 :     gk = cgetg(k, t_VEC); gel(gk,1) = gel(g,i);
    3265        4284 :     for (j=2; j<k; j++)
    3266        4214 :       gel(gk,j) = idealmoddivisor(bnr, idealmul(nf, gel(gk,j-1), gel(gk,1)));
    3267          70 :     gel(pows,i) = gk; /* powers of g[i] */
    3268             :   }
    3269             : 
    3270          70 :   C = cgetg(lc+1, t_VECSMALL); C[1] = c[lc];
    3271          70 :   for (i=2; i<=lc; i++) C[i] = C[i-1] * c[lc-i+1];
    3272             :   /* C[i] = c(k-i+1) * ... * ck */
    3273             :   /* j < C[i+1] <==> j only involves g(k-i)...gk */
    3274          70 :   i = 1;
    3275        4354 :   for (j=1; j < C[1]; j++)
    3276        4284 :     gel(list, j+1) = gmael(pows,lc,j);
    3277         140 :   while(j<no)
    3278             :   {
    3279             :     long k;
    3280             :     GEN a;
    3281           0 :     if (j == C[i+1]) i++;
    3282           0 :     a = gmael(pows,lc-i,j/C[i]);
    3283           0 :     k = j%C[i] + 1;
    3284           0 :     if (k > 1) a = idealmoddivisor(bnr, idealmul(nf, a, gel(list,k)));
    3285           0 :     gel(list, ++j) = a;
    3286             :   }
    3287          70 :   return list;
    3288             : }
    3289             : 
    3290             : /* x quadratic integer (approximate), recognize it. If error return NULL */
    3291             : static GEN
    3292        4424 : findbezk(GEN nf, GEN x)
    3293             : {
    3294        4424 :   GEN a,b, M = nf_get_M(nf), u = gcoeff(M,1,2);
    3295             :   long ea, eb;
    3296             : 
    3297             :   /* u t_COMPLEX generator of nf.zk, write x ~ a + b u, a,b in Z */
    3298        4424 :   b = grndtoi(mpdiv(imag_i(x), gel(u,2)), &eb);
    3299        4424 :   if (eb > -20) return NULL;
    3300        4424 :   a = grndtoi(mpsub(real_i(x), mpmul(b,gel(u,1))), &ea);
    3301        4424 :   if (ea > -20) return NULL;
    3302        4424 :   return signe(b)? coltoalg(nf, mkcol2(a,b)): a;
    3303             : }
    3304             : 
    3305             : static GEN
    3306          70 : findbezk_pol(GEN nf, GEN x)
    3307             : {
    3308          70 :   long i, lx = lg(x);
    3309          70 :   GEN y = cgetg(lx,t_POL);
    3310        4494 :   for (i=2; i<lx; i++)
    3311        4424 :     if (! (gel(y,i) = findbezk(nf,gel(x,i))) ) return NULL;
    3312          70 :   y[1] = x[1]; return y;
    3313             : }
    3314             : 
    3315             : /* P approximation computed at initial precision prec. Compute needed prec
    3316             :  * to know P with 1 word worth of trailing decimals */
    3317             : static long
    3318           0 : get_prec(GEN P, long prec)
    3319             : {
    3320           0 :   long k = gprecision(P);
    3321           0 :   if (k == 3) return precdbl(prec); /* approximation not trustworthy */
    3322           0 :   k = prec - k; /* lost precision when computing P */
    3323           0 :   if (k < 0) k = 0;
    3324           0 :   k += nbits2prec(gexpo(P) + 128);
    3325           0 :   if (k <= prec) k = precdbl(prec); /* dubious: old prec should have worked */
    3326           0 :   return k;
    3327             : }
    3328             : 
    3329             : /* Compute data for ellphist */
    3330             : static GEN
    3331        4354 : ellphistinit(GEN om, long prec)
    3332             : {
    3333        4354 :   GEN res,om1b,om2b, om1 = gel(om,1), om2 = gel(om,2);
    3334             : 
    3335        4354 :   if (gsigne(imag_i(gdiv(om1,om2))) < 0) { swap(om1,om2); om = mkvec2(om1,om2); }
    3336        4354 :   om1b = gconj(om1);
    3337        4354 :   om2b = gconj(om2); res = cgetg(4,t_VEC);
    3338        4354 :   gel(res,1) = gdivgs(elleisnum(om,2,0,prec),12);
    3339        4354 :   gel(res,2) = gdiv(PiI2(prec), gmul(om2, imag_i(gmul(om1b,om2))));
    3340        4354 :   gel(res,3) = om2b; return res;
    3341             : }
    3342             : 
    3343             : /* Computes log(phi^*(z,om)), using res computed by ellphistinit */
    3344             : static GEN
    3345        8708 : ellphist(GEN om, GEN res, GEN z, long prec)
    3346             : {
    3347        8708 :   GEN u = imag_i(gmul(z, gel(res,3)));
    3348        8708 :   GEN zst = gsub(gmul(u, gel(res,2)), gmul(z,gel(res,1)));
    3349        8708 :   return gsub(ellsigma(om,z,1,prec),gmul2n(gmul(z,zst),-1));
    3350             : }
    3351             : 
    3352             : /* Computes phi^*(la,om)/phi^*(1,om) where (1,om) is an oriented basis of the
    3353             :    ideal gf*gc^{-1} */
    3354             : static GEN
    3355        4354 : computeth2(GEN om, GEN la, long prec)
    3356             : {
    3357        4354 :   GEN p1,p2,res = ellphistinit(om,prec);
    3358             : 
    3359        4354 :   p1 = gsub(ellphist(om,res,la,prec), ellphist(om,res,gen_1,prec));
    3360        4354 :   p2 = imag_i(p1);
    3361        4354 :   if (gexpo(real_i(p1))>20 || gexpo(p2)> prec2nbits(minss(prec,realprec(p2)))-10)
    3362           0 :     return NULL;
    3363        4354 :   return gexp(p1,prec);
    3364             : }
    3365             : 
    3366             : /* Computes P_2(X)=polynomial in Z_K[X] closest to prod_gc(X-th2(gc)) where
    3367             :    the product is over the ray class group bnr.*/
    3368             : static GEN
    3369          70 : computeP2(GEN bnr, long prec)
    3370             : {
    3371          70 :   long clrayno, i, first = 1;
    3372          70 :   pari_sp av=avma, av2;
    3373          70 :   GEN listray, P0, P, lanum, la = get_lambda(bnr);
    3374          70 :   GEN nf = bnr_get_nf(bnr), f = gel(bnr_get_mod(bnr), 1);
    3375          70 :   listray = getallelts(bnr);
    3376          70 :   clrayno = lg(listray)-1; av2 = avma;
    3377             : PRECPB:
    3378          70 :   if (!first)
    3379             :   {
    3380           0 :     if (DEBUGLEVEL) pari_warn(warnprec,"computeP2",prec);
    3381           0 :     nf = gerepilecopy(av2, nfnewprec_shallow(checknf(bnr),prec));
    3382             :   }
    3383          70 :   first = 0; lanum = to_approx(nf,la);
    3384          70 :   P = cgetg(clrayno+1,t_VEC);
    3385        4424 :   for (i=1; i<=clrayno; i++)
    3386             :   {
    3387        4354 :     GEN om = get_om(nf, idealdiv(nf,f,gel(listray,i)));
    3388        4354 :     GEN s = computeth2(om,lanum,prec);
    3389        4354 :     if (!s) { prec = precdbl(prec); goto PRECPB; }
    3390        4354 :     gel(P,i) = s;
    3391             :   }
    3392          70 :   P0 = roots_to_pol(P, 0);
    3393          70 :   P = findbezk_pol(nf, P0);
    3394          70 :   if (!P) { prec = get_prec(P0, prec); goto PRECPB; }
    3395          70 :   return gerepilecopy(av, P);
    3396             : }
    3397             : 
    3398             : #define nexta(a) (a>0 ? -a : 1-a)
    3399             : static GEN
    3400          49 : do_compo(GEN A0, GEN B)
    3401             : {
    3402          49 :   long a, i, l = lg(B), v = fetch_var_higher();
    3403             :   GEN A, z;
    3404             :   /* now v > x = pol_x(0) > nf variable */
    3405          49 :   B = leafcopy(B); setvarn(B, v);
    3406          49 :   for (i = 2; i < l; i++) gel(B,i) = monomial(gel(B,i), l-i-1, 0);
    3407             :   /* B := x^deg(B) B(v/x) */
    3408          49 :   A = A0 = leafcopy(A0); setvarn(A0, v);
    3409          56 :   for  (a = 0;; a = nexta(a))
    3410             :   {
    3411          56 :     if (a) A = RgX_translate(A0, stoi(a));
    3412          56 :     z = resultant(A,B); /* in variable 0 */
    3413          56 :     if (issquarefree(z)) break;
    3414           7 :   }
    3415          49 :   (void)delete_var(); return z;
    3416             : }
    3417             : #undef nexta
    3418             : 
    3419             : static GEN
    3420          14 : galoisapplypol(GEN nf, GEN s, GEN x)
    3421             : {
    3422          14 :   long i, lx = lg(x);
    3423          14 :   GEN y = cgetg(lx,t_POL);
    3424             : 
    3425          14 :   for (i=2; i<lx; i++) gel(y,i) = galoisapply(nf,s,gel(x,i));
    3426          14 :   y[1] = x[1]; return y;
    3427             : }
    3428             : /* x quadratic, write it as ua + v, u,v rational */
    3429             : static GEN
    3430          70 : findquad(GEN a, GEN x, GEN p)
    3431             : {
    3432             :   long tu, tv;
    3433          70 :   pari_sp av = avma;
    3434             :   GEN u,v;
    3435          70 :   if (typ(x) == t_POLMOD) x = gel(x,2);
    3436          70 :   if (typ(a) == t_POLMOD) a = gel(a,2);
    3437          70 :   u = poldivrem(x, a, &v);
    3438          70 :   u = simplify_shallow(u); tu = typ(u);
    3439          70 :   v = simplify_shallow(v); tv = typ(v);
    3440          70 :   if (!is_scalar_t(tu)) pari_err_TYPE("findquad", u);
    3441          70 :   if (!is_scalar_t(tv)) pari_err_TYPE("findquad", v);
    3442          70 :   x = deg1pol(u, v, varn(a));
    3443          70 :   if (typ(x) == t_POL) x = gmodulo(x,p);
    3444          70 :   return gerepileupto(av, x);
    3445             : }
    3446             : static GEN
    3447          14 : findquad_pol(GEN p, GEN a, GEN x)
    3448             : {
    3449          14 :   long i, lx = lg(x);
    3450          14 :   GEN y = cgetg(lx,t_POL);
    3451          14 :   for (i=2; i<lx; i++) gel(y,i) = findquad(a, gel(x,i), p);
    3452          14 :   y[1] = x[1]; return y;
    3453             : }
    3454             : static GEN
    3455          49 : compocyclo(GEN nf, long m, long d)
    3456             : {
    3457          49 :   GEN sb,a,b,s,p1,p2,p3,res,polL,polLK,nfL, D = nf_get_disc(nf);
    3458             :   long ell,vx;
    3459             : 
    3460          49 :   p1 = quadhilbertimag(D);
    3461          49 :   p2 = polcyclo(m,0);
    3462          49 :   if (d==1) return do_compo(p1,p2);
    3463             : 
    3464          14 :   ell = m&1 ? m : (m>>2);
    3465          14 :   if (absequalui(ell,D)) /* ell = |D| */
    3466             :   {
    3467           0 :     p2 = gcoeff(nffactor(nf,p2),1,1);
    3468           0 :     return do_compo(p1,p2);
    3469             :   }
    3470          14 :   if (ell%4 == 3) ell = -ell;
    3471             :   /* nf = K = Q(a), L = K(b) quadratic extension = Q(t) */
    3472          14 :   polLK = quadpoly(stoi(ell)); /* relative polynomial */
    3473          14 :   res = rnfequation2(nf, polLK);
    3474          14 :   vx = nf_get_varn(nf);
    3475          14 :   polL = gsubst(gel(res,1),0,pol_x(vx)); /* = charpoly(t) */
    3476          14 :   a = gsubst(lift_shallow(gel(res,2)), 0,pol_x(vx));
    3477          14 :   b = gsub(pol_x(vx), gmul(gel(res,3), a));
    3478          14 :   nfL = nfinit(polL, DEFAULTPREC);
    3479          14 :   p1 = gcoeff(nffactor(nfL,p1),1,1);
    3480          14 :   p2 = gcoeff(nffactor(nfL,p2),1,1);
    3481          14 :   p3 = do_compo(p1,p2); /* relative equation over L */
    3482             :   /* compute non trivial s in Gal(L / K) */
    3483          14 :   sb= gneg(gadd(b, RgX_coeff(polLK,1))); /* s(b) = Tr(b) - b */
    3484          14 :   s = gadd(pol_x(vx), gsub(sb, b)); /* s(t) = t + s(b) - b */
    3485          14 :   p3 = gmul(p3, galoisapplypol(nfL, s, p3));
    3486          14 :   return findquad_pol(nf_get_pol(nf), a, p3);
    3487             : }
    3488             : 
    3489             : /* I integral ideal in HNF. (x) = I, x small in Z ? */
    3490             : static long
    3491         119 : isZ(GEN I)
    3492             : {
    3493         119 :   GEN x = gcoeff(I,1,1);
    3494         119 :   if (signe(gcoeff(I,1,2)) || !equalii(x, gcoeff(I,2,2))) return 0;
    3495         105 :   return is_bigint(x)? -1: itos(x);
    3496             : }
    3497             : 
    3498             : /* Treat special cases directly. return NULL if not special case */
    3499             : static GEN
    3500         119 : treatspecialsigma(GEN bnr)
    3501             : {
    3502         119 :   GEN bnf = bnr_get_bnf(bnr), nf = bnf_get_nf(bnf);
    3503         119 :   GEN f = gel(bnr_get_mod(bnr), 1),  D = nf_get_disc(nf);
    3504             :   GEN p1, p2;
    3505         119 :   long Ds, fl, tryf, i = isZ(f);
    3506             : 
    3507         119 :   if (i == 1) return quadhilbertimag(D); /* f = 1 */
    3508             : 
    3509         119 :   if (absequaliu(D,3)) /* Q(j) */
    3510             :   {
    3511           0 :     if (i == 4 || i == 5 || i == 7) return polcyclo(i,0);
    3512           0 :     if (!absequaliu(gcoeff(f,1,1),9) || !absequaliu(Q_content(f),3)) return NULL;
    3513             :     /* f = P_3^3 */
    3514           0 :     p1 = mkpolmod(bnf_get_tuU(bnf), nf_get_pol(nf));
    3515           0 :     return gadd(pol_xn(3,0), p1); /* x^3+j */
    3516             :   }
    3517         119 :   if (absequaliu(D,4)) /* Q(i) */
    3518             :   {
    3519          14 :     if (i == 3 || i == 5) return polcyclo(i,0);
    3520          14 :     if (i != 4) return NULL;
    3521           0 :     p1 = mkpolmod(bnf_get_tuU(bnf), nf_get_pol(nf));
    3522           0 :     return gadd(pol_xn(2,0), p1); /* x^2+i */
    3523             :   }
    3524         105 :   Ds = smodis(D,48);
    3525         105 :   if (i)
    3526             :   {
    3527          91 :     if (i==2 && Ds%16== 8) return compocyclo(nf, 4,1);
    3528          84 :     if (i==3 && Ds% 3== 1) return compocyclo(nf, 3,1);
    3529          70 :     if (i==4 && Ds% 8== 1) return compocyclo(nf, 4,1);
    3530          63 :     if (i==6 && Ds   ==40) return compocyclo(nf,12,1);
    3531          56 :     return NULL;
    3532             :   }
    3533             : 
    3534          14 :   p1 = gcoeff(f,1,1); /* integer > 0 */
    3535          14 :   tryf = itou_or_0(p1); if (!tryf) return NULL;
    3536          14 :   p2 = gcoeff(f,2,2); /* integer > 0 */
    3537          14 :   if (is_pm1(p2)) fl = 0;
    3538             :   else {
    3539           0 :     if (Ds % 16 != 8 || !absequaliu(Q_content(f),2)) return NULL;
    3540           0 :     fl = 1; tryf >>= 1;
    3541             :   }
    3542          14 :   if (tryf <= 3 || umodiu(D, tryf) || !uisprime(tryf)) return NULL;
    3543          14 :   if (fl) tryf <<= 2;
    3544          14 :   return compocyclo(nf,tryf,2);
    3545             : }
    3546             : 
    3547             : GEN
    3548         161 : quadray(GEN D, GEN f, long prec)
    3549             : {
    3550             :   GEN bnr, y, bnf;
    3551         161 :   pari_sp av = avma;
    3552             : 
    3553         161 :   if (isint1(f)) return quadhilbert(D, prec);
    3554         126 :   quadray_init(&D, f, &bnf, prec);
    3555         126 :   bnr = Buchray(bnf, f, nf_INIT|nf_GEN);
    3556         126 :   if (is_pm1(bnr_get_no(bnr))) { avma = av; return pol_x(0); }
    3557         126 :   if (signe(D) > 0)
    3558           7 :     y = bnrstark(bnr,NULL,prec);
    3559             :   else
    3560             :   {
    3561         119 :     bnr = gel(bnrconductor_i(bnr,NULL,2), 2);
    3562         119 :     y = treatspecialsigma(bnr);
    3563         119 :     if (!y) y = computeP2(bnr, prec);
    3564             :   }
    3565         126 :   return gerepileupto(av, y);
    3566             : }

Generated by: LCOV version 1.11