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

Generated by: LCOV version 1.16