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 - basemath - grossenchar.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.14.0 lcov report (development 27775-aca467eab2) Lines: 1032 1046 98.7 %
Date: 2022-07-03 07:33:15 Functions: 65 65 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2000  The PARI group.
       2             : 
       3             :    This file is part of the PARI/GP package.
       4             : 
       5             :    PARI/GP is free software; you can redistribute it and/or modify it under the
       6             :    terms of the GNU General Public License as published by the Free Software
       7             :    Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             :    ANY WARRANTY WHATSOEVER.
       9             : 
      10             :    Check the License for details. You should have received a copy of it, along
      11             :    with the package; see the file 'COPYING'. If not, write to the Free Software
      12             :    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : #include "pari.h"
      14             : #include "paripriv.h"
      15             : 
      16             : #define DEBUGLEVEL DEBUGLEVEL_gchar
      17             : 
      18             : static GEN gchari_eval(GEN gc, GEN chi, GEN x, long flag, GEN logchi, GEN s0, long prec);
      19             : 
      20             : /*********************************************************************/
      21             : /**                                                                 **/
      22             : /**                    HECKE GROSSENCHARACTERS                      **/
      23             : /**         contributed by Pascal Molin and Aurel Page (2022)       **/
      24             : /**                                                                 **/
      25             : /*********************************************************************/
      26             : 
      27             : /*
      28             :   characters can be represented by:
      29             :    - t_COL of coordinates on the snf basis (mostly for gp use): prefix gchar_
      30             :    - t_VEC of coordinates on the internal basis: prefix gchari_
      31             :    - t_VEC of R/Z components (logs): prefix gcharlog_
      32             : 
      33             :    see gchar_internal for SNF -> internal
      34             :    and gchari_duallog for internal -> R/Z components
      35             : */
      36             : 
      37             : /*
      38             : localstar: represent (Z_F/m)^* . {+-1}^r + generators of U_{i-1}(p)/U_i
      39             : structure:
      40             : - [ sprk for p^k | m ] , size np
      41             : - [ Ufil_p for p | m ], size np
      42             : - m_oo, t_VECSMALL of size nci <= r1 (indices of real places)
      43             : 
      44             : where Ufil_p = [ Mat([gen[j], t_COL of size ncp]_j) ]_{1<=i<=k}
      45             : */
      46             : 
      47             : #define GC_LENGTH   12
      48             : #define LOCS_LENGTH 4
      49             : 
      50             : static GEN
      51         532 : compute_Lcyc(GEN Lsprk, GEN moo)
      52             : {
      53         532 :   long i, l = lg(Lsprk), len = l+lg(moo)-1;
      54         532 :   GEN Lcyc = cgetg(len,t_VEC);
      55        1043 :   for (i = 1; i < l; i++)   gel(Lcyc,i) = sprk_get_cyc(gel(Lsprk,i));
      56         651 :   for (     ; i < len; i++) gel(Lcyc,i) = mkvec(gen_2);
      57         532 :   return Lcyc;
      58             : }
      59             : 
      60             : static long
      61        1757 : sprk_get_ncp(GEN sprk) { return lg(sprk_get_cyc(sprk)) - 1; }
      62             : 
      63             : /* true nf; modulus = [ factor(m_f), m_oo ] */
      64             : static GEN
      65         532 : localstar(GEN nf, GEN famod, GEN moo)
      66             : {
      67         532 :   GEN Lcyc, cyc, Lsprk, Lgenfil, P = gel(famod,1), E = gel(famod,2);
      68         532 :   long i, l = lg(P);
      69             : 
      70         532 :   Lsprk = cgetg(l, t_VEC);
      71         532 :   Lgenfil = cgetg(l, t_VEC);
      72        1043 :   for (i = 1; i < l; i++)
      73             :   {
      74         511 :     long e, k = itos(gel(E,i));
      75         511 :     GEN v, sprk = log_prk_init(nf, gel(P,i), k, NULL);
      76         511 :     gel(Lsprk,i) = sprk;
      77         511 :     gel(Lgenfil,i) = v = cgetg(k+1, t_VEC);
      78             :     /* log on sprk of generators of U_{e-1}/U_e(pr) */
      79         511 :     gel(v, 1) = col_ei(sprk_get_ncp(sprk), 1);
      80        1113 :     for (e = 2; e <= k; e++) gel(v, e) = sprk_log_gen_pr2(nf, sprk, e);
      81             :   }
      82         532 :   Lcyc = compute_Lcyc(Lsprk, moo);
      83         532 :   if (lg(Lcyc) > 1)
      84         252 :     cyc = shallowconcat1(Lcyc);
      85             :   else
      86         280 :     cyc = cgetg(1, t_VEC);
      87         532 :   return mkvec4(cyc, Lsprk, Lgenfil, mkvec2(famod,moo));
      88             : }
      89             : 
      90             : /* (nv * log|x^sigma|/norm, arg(x^sigma))/2*Pi
      91             :  * substract norm so that we project to the hyperplane
      92             :  * H : sum n_s x_s = 0 */
      93             : static GEN
      94      194074 : nfembedlog(GEN *pnf, GEN x, long prec)
      95             : {
      96      194074 :   pari_sp av = avma;
      97      194074 :   GEN logs, cxlogs, nf = *pnf;
      98      194074 :   long k, r1, r2, n, extrabit, extranfbit = 0, nfprec, nfprec0, logprec;
      99             : 
     100      194074 :   nfprec0 = nf_get_prec(nf);
     101      194074 :   nf_get_sign(nf, &r1, &r2);
     102      194074 :   n = r1 + 2*r2;
     103      194074 :   logprec = prec;
     104      194074 :   extrabit = expu(n) + gexpo(nf_get_M(nf)) + 10;
     105      194074 :   if (typ(x) == t_MAT)
     106             :   {
     107      193542 :     long l = lg(gel(x,2));
     108      193542 :     if (l > 1)
     109             :     {
     110      193542 :       extrabit += expu(l) + gexpo(gel(x,2));
     111      193542 :       extranfbit = gexpo(gel(x,1));
     112             :     }
     113             :   }
     114             :   else
     115             :   {
     116         532 :     x = nf_to_scalar_or_basis(nf,x);
     117         532 :     extranfbit = gexpo(x);
     118             :   }
     119      194074 :   if (DEBUGLEVEL>3)
     120           0 :     err_printf("  nfembedlog: prec=%d extrabit=%d nfprec=%d extralogprec=%d\n",
     121             :                prec, nbits2extraprec(extrabit + extranfbit), nfprec0,
     122             :                nbits2extraprec(extrabit));
     123      194074 :   nfprec = prec + nbits2extraprec(extrabit + extranfbit);
     124      194074 :   logprec = prec + nbits2extraprec(extrabit);
     125      194074 :   if (nfprec0 < nfprec)
     126             :   {
     127        1047 :     if (DEBUGLEVEL)
     128           0 :       err_printf("  nfembedlog: increasing prec %d -> %d\n", nfprec0, nfprec);
     129        1047 :     *pnf = nf = nfnewprec_shallow(nf, nfprec);
     130        1047 :     av = avma;
     131             :   }
     132      194074 :   if (!(cxlogs = nf_cxlog(nf, x, logprec))) return gc_NULL(av);
     133      194060 :   if (!(cxlogs = nf_cxlog_normalize(nf, cxlogs, logprec))) return gc_NULL(av);
     134      194060 :   logs = cgetg(n+1,t_COL);
     135      717986 :   for (k = 1; k <= r1+r2; k++) gel(logs,k) = real_i(gel(cxlogs,k));
     136      408983 :   for (     ; k <= n; k++) gel(logs,k) = gmul2n(imag_i(gel(cxlogs,k-r2)), -1);
     137      194060 :   extrabit = gexpo(logs);
     138      194060 :   if (extrabit < 0) extrabit = 0;
     139      194060 :   prec += nbits2extraprec(extrabit);
     140      194060 :   return gerepileupto(av, gdiv(logs, Pi2n(1,prec)));
     141             : }
     142             : 
     143             : static GEN
     144        1736 : gchar_Sval(GEN nf, GEN S, GEN x)
     145             : {
     146        1736 :   GEN res = cgetg(lg(S),t_COL);
     147             :   long i;
     148        1736 :   if (typ(x)==t_MAT)
     149        3381 :     for (i=1; i<lg(S); i++)
     150        2177 :       gel(res, i) = famat_nfvalrem(nf, x, gel(S,i), NULL);
     151             :   else
     152         896 :     for (i=1; i<lg(S); i++)
     153         364 :       gel(res, i) = stoi(nfval(nf, x, gel(S,i)));
     154        1736 :   return res;
     155             : }
     156             : 
     157             : /* true nf; log_prk(x*pi_pr^{-v_pr(x)}), sign(sigma(x)) */
     158             : static GEN
     159      189581 : gchar_logm(GEN nf, GEN locs, GEN x)
     160             : {
     161      189581 :   GEN moo, loga, Lsprk = locs_get_Lsprk(locs);
     162      189581 :   long i, l = lg(Lsprk);
     163             : 
     164      189581 :   moo = locs_get_m_infty(locs);
     165      189581 :   if (typ(x) != t_MAT) x = to_famat_shallow(x, gen_1);
     166      189581 :   loga = cgetg(l+1, t_VEC);
     167      606522 :   for (i = 1; i < l; i++)
     168             :   {
     169      416941 :     GEN sprk = gel(Lsprk, i), pr = sprk_get_pr(sprk);
     170      416941 :     GEN g = vec_append(gel(x,1), pr_get_gen(pr));
     171      416941 :     GEN v = famat_nfvalrem(nf, x, pr, NULL);
     172      416941 :     GEN e = vec_append(gel(x,2), gneg(v));
     173      416941 :     gel(loga, i) = famat_zlog_pr(nf, g, e, sprk, NULL);
     174             :   }
     175      189581 :   gel(loga, i) = zc_to_ZC( nfsign_arch(nf, x, moo) );
     176      189581 :   return shallowconcat1(loga);
     177             : }
     178             : 
     179             : static GEN
     180        1750 : gchar_nflog(GEN *pnf, GEN zm, GEN S, GEN x, long prec)
     181             : {
     182        1750 :   GEN emb = nfembedlog(pnf, x, prec);
     183        1750 :   if (!emb) return NULL;
     184        1736 :   return shallowconcat1(mkvec3(gchar_Sval(*pnf,S,x),
     185             :                                gchar_logm(*pnf,zm,x), emb));
     186             : }
     187             : 
     188             : /*******************************************************************/
     189             : /*                                                                 */
     190             : /*                        CHARACTER GROUP                          */
     191             : /*                                                                 */
     192             : /*******************************************************************/
     193             : 
     194             : /* embed v in vector of length size, shifted to the right */
     195             : static GEN
     196        6097 : embedcol(GEN v, long size, long shift)
     197             : {
     198             :   long k;
     199        6097 :   GEN w = zerocol(size);
     200       48930 :   for (k = 1; k < lg(v); k++) gel(w, shift+k) = gel(v,k);
     201        6097 :   return w;
     202             : }
     203             : 
     204             : /* write exact rationals from real approximation, in place. */
     205             : static void
     206        4658 : shallow_clean_rat(GEN v, long k0, long k1, GEN den, long prec)
     207             : {
     208        4658 :   long k, e, bit = -prec2nbits(prec)/2;
     209       43438 :   for (k = k0; k <= k1; k++)
     210             :   {
     211       38780 :     GEN t = gel(v,k);
     212       38780 :     if (den) t = gmul(t, den);
     213       38780 :     t = grndtoi(t, &e);
     214       38780 :     if (DEBUGLEVEL>1) err_printf("[%Ps*%Ps=%Ps..e=%ld|prec=%ld]\n",
     215           0 :                                  gel(v,k), den? den: gen_1, t, e, prec);
     216       38780 :     if (e > bit)
     217             :       pari_err_BUG("gcharinit, non rational entry"); /*LCOV_EXCL_LINE*/
     218       38780 :     gel(v, k) = den? gdiv(t, den): t;
     219             :   }
     220        4658 : }
     221             : 
     222             : /* FIXME: export ? */
     223             : static GEN
     224        1050 : rowreverse(GEN m)
     225             : {
     226             :   long i, l;
     227             :   GEN v;
     228        1050 :   if (lg(m) == 1) return m;
     229        1050 :   l = lgcols(m); v = cgetg(l, t_VECSMALL);
     230        4025 :   for (i = 1; i < l; i++) v[i] = l - i;
     231        1050 :   return rowpermute(m, v);
     232             : }
     233             : 
     234             : /* lower-left hnf on subblock m[r0+1..r0+nr, c0+1..c0+nc]
     235             :  * return base change matrix U */
     236             : static GEN
     237        1050 : hnf_block(GEN m, long r0, long nr, long c0, long nc)
     238             : {
     239             :   GEN block, u, uu;
     240        1050 :   long nm = lg(m)-1, k;
     241        1050 :   pari_sp av = avma;
     242             : 
     243        1050 :   block = matslice(m, r0+1, r0+nr, c0+1, c0+nc);
     244        1050 :   block = Q_remove_denom(block, NULL);
     245        1050 :   block = rowreverse(block);
     246             : 
     247        1050 :   (void)ZM_hnfall(block, &u, 1);
     248        1050 :   vecreverse_inplace(u);
     249        1050 :   uu = matid(nm);
     250             :   /* embed in matid */
     251        5502 :   for (k = 1; k <= nc; k++)
     252        4452 :     gel(uu,c0+k) = embedcol(gel(u,k),nm,c0);
     253        1050 :   return gerepilecopy(av, uu);
     254             : }
     255             : 
     256             : static GEN
     257         647 : lll_block(GEN m, long r0, long nr, long c0, long nc)
     258             : {
     259             :   GEN block, u, uu;
     260         647 :   long nm = lg(m)-1, k;
     261         647 :   pari_sp av = avma;
     262             : 
     263         647 :   block = matslice(m, r0+1, r0+nr, c0+1, c0+nc);
     264         647 :   u = lll(block); vecreverse_inplace(u);
     265         647 :   if (lg(u) <= nc) return NULL;
     266         644 :   uu = matid(nm); /* embed in matid */
     267        2289 :   for (k = 1; k <= nc; k++) gel(uu,c0+k) = embedcol(gel(u,k),nm,c0);
     268         644 :   return gerepilecopy(av, uu);
     269             : }
     270             : 
     271             : /* insert a column at index i, no copy */
     272             : static GEN
     273        2329 : shallowmatinsert(GEN m, GEN x, long i)
     274             : {
     275        2329 :   long k, n = lg(m);
     276        2329 :   GEN mm = cgetg(n+1,t_MAT);
     277       14235 :   for (k=1; k < i; k++) gel(mm, k) = gel(m, k);
     278        2329 :   gel(mm, i) = x;
     279        3902 :   for (k=i; k < n; k++) gel(mm, k+1) = gel(m, k);
     280        2329 :   return mm;
     281             : }
     282             : 
     283             : static GEN
     284        2329 : vec_v0(long n, long n0, long r1, long r2)
     285             : {
     286             :   long k;
     287        2329 :   GEN C = zerocol(n);
     288        5394 :   for (k = 1; k <= r1; k++) gel(C, n0++) = gen_1;
     289        5641 :   for (k = 1; k <= r2; k++) gel(C, n0++) = gen_2;
     290        2329 :   return C;
     291             : }
     292             : 
     293             : /* select cm embeddings; return a matrix */
     294             : /* TODO detect if precision was insufficient */
     295             : static GEN
     296         203 : cm_select(GEN nf, GEN cm, long prec)
     297             : {
     298             :   GEN emb, keys, v, m_sel, imag_emb;
     299         203 :   long nc, d_cm, r_cm, c, i, j, r2 = nf_get_r2(nf);
     300             :   pari_sp av;
     301             : 
     302         203 :   d_cm = degpol(gel(cm, 1)); /* degree of the cm field; even */
     303         203 :   nc = d_cm / 2; /* nb of clusters */
     304         203 :   r_cm = nf_get_degree(nf) / d_cm; /* nb by cluster; nc * r_cm = r2 */
     305         203 :   m_sel = zeromatcopy(nc, r2); /* selection matrix */
     306         203 :   av = avma;
     307             :   /* group complex embeddings */
     308         203 :   emb = nfeltembed(nf, gel(cm, 2), NULL, prec);
     309             :   /* sort */
     310         203 :   imag_emb = imag_i(emb);
     311         203 :   keys = gadd(gmul(mppi(prec), real_i(emb)), gabs(imag_emb, prec));
     312         203 :   v = indexsort(keys);
     313             : 
     314         567 :   for (j = c = 1; c <= nc; c++)
     315             :   {
     316         364 :     int ref = gsigne(gel(imag_emb, v[j]));
     317         364 :     gcoeff(m_sel, c, v[j]) = gen_1;
     318         364 :     j++;
     319         476 :     for (i = 2; i <= r_cm; i++)
     320             :     {
     321         112 :       int s = gsigne(gel(imag_emb, v[j]));
     322         112 :       gcoeff(m_sel, c, v[j]) = (s == ref) ? gen_1 : gen_m1;
     323         112 :       j++;
     324             :     }
     325             :   }
     326         203 :   return gc_const(av, m_sel);
     327             : }
     328             : 
     329             : static GEN gchar_hnfreduce_shallow(GEN gc, GEN cm);
     330             : static void gchar_snfbasis_shallow(GEN gc, GEN rel);
     331             : static void gcharmat_tinverse(GEN gc, GEN m, long prec);
     332             : static GEN gcharmatnewprec_shallow(GEN gc, long mprec);
     333             : 
     334             : /* return a set S of prime ideals such that Cl_S(K) = 1 */
     335             : static GEN
     336         532 : bestS(GEN bnf)
     337             : {
     338         532 :   GEN v, S, hw, hv = bnf_get_no(bnf), DL, dl;
     339             :   long i, lS;
     340             :   ulong l;
     341             :   forprime_t P;
     342             : 
     343         532 :   if (equali1(hv)) return mkvec2(cgetg(1,t_VEC), cgetg(1,t_VEC));
     344         203 :   v = diagonal_shallow(bnf_get_cyc(bnf));
     345         203 :   S = cgetg(expi(hv)+2, t_VEC); lS = 1;
     346         203 :   DL = cgetg(expi(hv)+2, t_VEC);
     347         203 :   u_forprime_init(&P,2,ULONG_MAX);
     348         658 :   while ((l = u_forprime_next(&P)))
     349             :   {
     350         658 :     pari_sp av = avma;
     351         658 :     GEN w, Sl = idealprimedec(bnf, utoi(l));
     352         658 :     long nSl = lg(Sl)-1;
     353        1148 :     for (i = 1; i < nSl; i++) /* remove one prime ideal */
     354             :     {
     355         693 :       dl = isprincipal(bnf, gel(Sl,i));
     356         693 :       w = ZM_hnf(shallowconcat(v, dl));
     357         693 :       hw = ZM_det(w);
     358         693 :       if (cmpii(hw, hv) < 0)
     359             :       {
     360         364 :         gel(DL,lS) = dl;
     361         364 :         gel(S,lS++) = gel(Sl,i);
     362         364 :         hv = hw; v = w; av = avma;
     363         364 :         if (equali1(hv)) { setlg(S, lS); setlg(DL, lS); return mkvec2(S,DL); }
     364             :       }
     365             :     }
     366         455 :     set_avma(av);
     367             :   }
     368             :   return NULL;/*LCOV_EXCL_LINE*/
     369             : }
     370             : 
     371             : static GEN
     372         532 : gcharDLdata(GEN bnf, GEN S, GEN DL)
     373             : {
     374         532 :   GEN M, h, Minv, Lalpha, t, dl, alpha, gen, cyc = bnf_get_cyc(bnf);
     375             :   long i;
     376         532 :   M = shallowmatconcat(DL);
     377         532 :   h = bnf_get_no(bnf);
     378         532 :   gen = bnf_get_gen(bnf);
     379             : 
     380             :   /* compute right inverse of M modulo cyc */
     381         532 :   M = shallowtrans(M);
     382         532 :   M = shallowmatconcat(mkcol2(M,diagonal_shallow(cyc)));
     383         532 :   Minv = matinvmod(M,h);
     384         532 :   Minv = vecslice(Minv,1,lg(Minv)-lg(cyc));
     385         532 :   Minv = shallowtrans(Minv);
     386             : 
     387         532 :   Lalpha = cgetg(lg(Minv),t_VEC);
     388         847 :   for (i=1; i<lg(Minv); i++)
     389             :   {
     390             :     /* gen[i] = (alpha) * prod_j S[j]^Minv[j,i] */
     391         315 :     t = isprincipalfact(bnf, gel(gen,i), S, gneg(gel(Minv,i)), nf_GENMAT);
     392         315 :     dl = gel(t, 1); alpha = gel(t, 2);
     393         315 :     if (!gequal0(dl)) pari_err_BUG("gcharDLdata (non-principal ideal)");
     394         315 :     gel(Lalpha,i) = alpha;
     395             :   }
     396         532 :   return mkvec2(Minv, Lalpha);
     397             : }
     398             : 
     399             : /* compute basis of characters; gc[1] generating family, as rows */
     400             : GEN
     401         553 : gcharinit(GEN bnf, GEN mod, long prec)
     402             : {
     403         553 :   pari_sp av = avma;
     404             :   GEN nf, zm, zmcyc, S, DLdata, sfu, logx;
     405             :   GEN fa2, archp, z, C, gc, cm, cyc, rel, U, Ui, m, m_inv, m0, u0;
     406             :   long n, k, r1, r2, ns, nc, nu, nm, order;
     407         553 :   long evalprec = prec, nfprec, mprec, embprec;
     408             : 
     409         553 :   prec = evalprec + 1; /* default 1 extra word */
     410             : 
     411             :   /* note on precision:
     412             : 
     413             :      - evalprec = precision requested for evaluation
     414             : 
     415             :      - prec = precision available
     416             :             = (relative) precision of parameters = m_inv
     417             :             = (relative) precision of evaluation of small chars
     418             :               if no cancellation
     419             : 
     420             :      - nfprec = internal nf precision used for
     421             :        the embedding matrix m
     422             : 
     423             :      In the structure we store [evalprec,prec,nfprec]
     424             : 
     425             :      When evaluating chi(x) at evalprec,
     426             :      we need prec >= max(evalprec + exponent(chi), nfprec+exponent(x))
     427             :      where exponent(x) is the exponent of the number field element alpha
     428             :      obtained after principalisation of x.
     429             : 
     430             :      If prec is not sufficient, we call gcharnewprec.
     431             : 
     432             :      To mitigate precision increase, we always initialize the structure
     433             :      with one extra word.
     434             : 
     435             :      Final remark: a gchar struct may have inconsistent values
     436             :      for prec and nfprec, for example if it has been saved and loaded at
     437             :      default prec : one should call gcharnewprec then.
     438             :   */
     439             : 
     440         553 :   if (!checkbnf_i(bnf))
     441             :   {
     442          84 :     nfprec = prec;
     443          84 :     bnf = bnfinit0(bnf, 1, NULL, nfprec);
     444          77 :     nf = shallowcopy(bnf_get_nf(bnf));
     445             :   }
     446             :   else
     447             :   {
     448         469 :     GEN fu = bnf_get_sunits(bnf);
     449         469 :     if (!fu) fu = bnf_get_fu(bnf); /* impose fundamental units */
     450         469 :     nf = shallowcopy(bnf_get_nf(bnf));
     451         469 :     nfprec = nf_get_prec(nf);
     452             :   }
     453             : 
     454             :   /* Dirichlet group + make sure mod contains archimedean places */
     455         546 :   mod = check_mod_factored(nf,mod,NULL,&fa2,&archp,NULL);
     456         532 :   sort_factor(fa2, (void*)&cmp_prime_ideal, &cmp_nodata);
     457         532 :   zm = localstar(nf, fa2, archp);
     458         532 :   zmcyc = locs_get_cyc(zm);
     459             : 
     460             :   /* set of primes S and valuations of generators */
     461         532 :   S = bestS(bnf);
     462         532 :   DLdata = gel(S,2);
     463         532 :   S = gel(S,1);
     464         532 :   DLdata = gcharDLdata(bnf, S, DLdata);
     465             : 
     466         532 :   nf_get_sign(nf, &r1, &r2);
     467         532 :   n = r1+2*r2;
     468         532 :   ns = lg(S) - 1;
     469         532 :   nu = r1+r2-1 + ns;
     470         532 :   nc = lg(zmcyc) - 1;
     471         532 :   nm = ns+nc+n; /* number of parameters = ns + nc + r1 + r2 + r2 */
     472             : 
     473             :   /* units and S-units */
     474         532 :   sfu = gel(bnfunits(bnf,S), 1);
     475         532 :   sfu = vec_shorten(sfu, nu); /* remove torsion */
     476             : 
     477             :   /* root of unity */
     478         532 :   order = bnf_get_tuN(bnf);
     479         532 :   z = bnf_get_tuU(bnf);
     480             : 
     481             :   /* maximal cm subfield */
     482         532 :   cm = nfsubfieldscm(nf, 0);
     483             : 
     484             :   /*
     485             :    Now compute matrix of parameters,
     486             :    until we obtain the right precision
     487             :    FIXME: make sure generators, units, and discrete logs
     488             :           do not depend on precision.
     489             : 
     490             :    m0 is the matrix of units embeddings
     491             :    u  is the HNF base change, m = m0*u
     492             : 
     493             :    subsequent steps may lead to precision increase, we put everything in gc
     494             :    struct and modify it in place.
     495             : 
     496             :      A) sets m0
     497             :      B) sets U, cyc, rel, U and Ui
     498             :      C) sets m_inv
     499             :   */
     500             : 
     501             :   /* A) make big matrix m0 of embeddings of units */
     502             : 
     503         532 :   if (DEBUGLEVEL>2) err_printf("start matrix m\n");
     504         532 :   m = cgetg(nm + 1, t_MAT);
     505         532 :   mprec = 3 + nbits2extraprec(nm+10);
     506         532 :   embprec = mprec;
     507             :   for(;;)
     508             :   {
     509        1750 :     for (k = 1; k <= nu; k++)
     510             :     { /* Lambda_S (S-units) then Lambda_f, fund. units */
     511        1218 :       logx = gchar_nflog(&nf, zm, S, gel(sfu,k), embprec);
     512        1218 :       if (!logx) break;
     513        1204 :       gel(m, k) = logx;
     514             :     }
     515         546 :     if (k > nu) break;
     516          14 :     if (DEBUGLEVEL) err_printf("gcharinit: increasing embprec %d -> %d\n",
     517             :                                embprec, precdbl(embprec));
     518          14 :     embprec = precdbl(embprec);
     519             :   }
     520        1435 :   for (k = 1; k <= nc; k++) /* Gamma, structure of (Z/m)* */
     521             :   {
     522         903 :     C = zerocol(nm);
     523         903 :     gel(C, ns+k) = gel(zmcyc, k);
     524         903 :     gel(m, nu+k) = C;
     525             :   }
     526             :   /* zeta, root of unity */
     527         532 :   gel(m, nu+nc+1) = gchar_nflog(&nf, zm, S, z, mprec);
     528         532 :   shallow_clean_rat(gel(m, nu+nc+1), 1, nm, stoi(order), mprec);
     529        1337 :   for (k = 1; k <= r2; k++) /* embed Z^r_2 */
     530             :   {
     531         805 :     C = zerocol(nm);
     532         805 :     gel(C, ns+nc+r1+r2+k) = gen_1;
     533         805 :     gel(m, nu+nc+1+k) = C;
     534             :   }
     535         532 :   if (DEBUGLEVEL>1) err_printf("matrix m = %Ps\n", m);
     536             : 
     537         532 :   m0 = m;
     538         532 :   u0 = gen_0;
     539         532 :   rel = U = Ui = gen_0;
     540         532 :   cyc = gen_0;
     541         532 :   m_inv = gen_0;
     542             : 
     543             :   /* only m and m_inv depend on prec and are recomputed under gcharnewprec. */
     544         532 :   gc = mkvecn(GC_LENGTH,
     545             :               m_inv, /* internal basis, characters as rows */
     546             :               bnf,
     547             :               nf,
     548             :               zm,    /* Zk/mod, nc components */
     549             :               S,     /* generators of clgp, ns components */
     550             :               DLdata,
     551             :               sfu,
     552             :               mkvec2(mkvecsmall3(evalprec,prec,nfprec),
     553             :                      mkvecsmall3(0,0,0)), /* ntors, nfree, nalg */
     554             :               cyc, /* reduced components */
     555             :               mkvec3(rel, U, Ui), /* internal / SNF base change */
     556             :               m0,                 /* embeddings of units */
     557             :               u0);                /* m_inv = (m0 u0)~^-1 */
     558             : 
     559             :   /* B) do HNF reductions + LLL (may increase precision) */
     560         532 :   m = gchar_hnfreduce_shallow(gc, cm);
     561             : 
     562             :   /* C) compute snf basis of torsion subgroup */
     563         532 :   rel = shallowtrans(matslice(m, 1, ns+nc, 1, ns+nc));
     564         532 :   gchar_snfbasis_shallow(gc, rel);
     565             : 
     566             :   /* D) transpose inverse m_inv = (m0*u)~^-1 (may increase precision) */
     567         532 :   gcharmat_tinverse(gc, m, prec);
     568         532 :   return gerepilecopy(av, gc);
     569             : }
     570             : 
     571             : /* b) do HNF reductions + LLL, keep base change u0 */
     572             : static GEN
     573         532 : gchar_hnfreduce_shallow(GEN gc, GEN cm)
     574             : {
     575         532 :   GEN bnf = gchar_get_bnf(gc), nf = gchar_get_nf(gc), u, u0, m;
     576         532 :   long order, r1, r2, ns, nc, n, nu, nm, ncm = 0, mprec;
     577             : 
     578         532 :   nf_get_sign(nf, &r1, &r2);
     579         532 :   n = r1 + 2*r2;
     580         532 :   nu = r1 + r2 - 1;
     581         532 :   ns = gchar_get_ns(gc);
     582         532 :   nc = gchar_get_nc(gc);
     583         532 :   nm = ns+nc+n; /* ns + nc + r1 + r2 + r2 */
     584         532 :   order = 2*bnf_get_tuN(bnf);
     585         532 :   u0 = matid(nm);
     586         532 :   m = shallowcopy(gchar_get_m0(gc)); /* keep m0 unchanged */
     587         532 :   mprec = gprecision(m);
     588         532 :   if (DEBUGLEVEL>1) err_printf("matrix m = %Ps\n", m);
     589         532 :   if (nc)
     590             :   { /* keep steps 1&2 to make sure we have zeta_m */
     591         252 :     u = hnf_block(m, ns,nc, ns+nu,nc+1);
     592         252 :     u0 = ZM_mul(u0, u); m = RgM_ZM_mul(m, u);
     593         252 :     if (DEBUGLEVEL>2) err_printf("step 1 -> %Ps\n", m);
     594         252 :     u = hnf_block(m, ns,nc, ns,nu+nc);
     595         252 :     u0 = ZM_mul(u0, u); m = RgM_ZM_mul(m, u);
     596         252 :     if (DEBUGLEVEL>2) err_printf("step 2 -> %Ps\n", m);
     597             :   }
     598         532 :   if (r2)
     599             :   {
     600         343 :     u = hnf_block(m, nm-r2,r2, nm-r2-1,r2+1);
     601         343 :     u0 = ZM_mul(u0, u); m = RgM_ZM_mul(m, u);
     602         343 :     if (DEBUGLEVEL>2) err_printf("step 3 -> %Ps\n", m);
     603             :   }
     604             :   /* remove last column */
     605         532 :   setlg(u0, nm); setlg(m, nm);
     606         532 :   if (DEBUGLEVEL>2) err_printf("remove last col -> %Ps\n", m);
     607             : 
     608         532 :   if (!gequal0(cm))
     609             :   {
     610             :     GEN v, Nargs;
     611         203 :     long bit = - prec2nbits(mprec) + 16 + expu(order);
     612             :     /* reduce on Norm arguments */
     613         203 :     v = cm_select(nf, cm, gchar_get_nfprec(gc));
     614         203 :     if (DEBUGLEVEL>2) err_printf("cm_select -> %Ps\n", v);
     615         203 :     ncm = nbrows(v);
     616         203 :     gchar_set_u0(gc, u0);
     617             :     for(;;)
     618          21 :     {
     619             :       long e, emax, i;
     620         224 :       Nargs = gmul(v, rowslice(m, nm-r2+1, nm));
     621         224 :       if (DEBUGLEVEL>2) err_printf("Nargs -> %Ps\n", Nargs);
     622         224 :       emax = bit-1;
     623        1094 :       for (i = ns+nc+1; i < lg(Nargs); i++)
     624             :       {
     625         870 :         gel(Nargs,i) = grndtoi(gmulgs(gel(Nargs,i), order), &e);
     626         870 :         emax = maxss(emax,e);
     627             :       }
     628         224 :       if (emax < bit) break;
     629          21 :       if (DEBUGLEVEL>1) err_printf("cm select: doubling prec\n");
     630          21 :       mprec = precdbl(mprec);
     631          21 :       m = gcharmatnewprec_shallow(gc, mprec);
     632             :     }
     633         203 :     if (DEBUGLEVEL>2) err_printf("rounded Nargs -> %Ps\n", Nargs);
     634         203 :     u = hnf_block(Nargs, 0, ncm, ns+nc, n-1);
     635         203 :     u0 = ZM_mul(u0, u); m = RgM_ZM_mul(m, u);
     636         203 :     if (DEBUGLEVEL>2) err_printf("after cm reduction -> %Ps\n", m);
     637             :   }
     638             : 
     639             :   /* apply LLL on Lambda_m, may need to increase prec */
     640         532 :   gchar_set_nalg(gc, ncm);
     641         532 :   gchar_set_u0(gc, u0);
     642             : 
     643             :   /* TODO factor these two LLL reduction codes in a function? */
     644         532 :   if (ncm > 0)
     645             :   {
     646         203 :     GEN u = NULL;
     647             :     while (1)
     648             :     {
     649         203 :       u = lll_block(m, ns+nc, n, ns+nc, ncm); if (u) break;
     650           0 :       mprec = precdbl(mprec);
     651             :       /* recompute m0 * u0 to higher prec */
     652           0 :       m = gcharmatnewprec_shallow(gc, mprec);
     653             :     }
     654         203 :     u0 = ZM_mul(u0, u); m = RgM_ZM_mul(m, u);
     655         203 :     if (DEBUGLEVEL>1) err_printf("after LLL reduction (CM block) -> %Ps\n", m);
     656             :   }
     657         532 :   gchar_set_u0(gc, u0);
     658             : 
     659         532 :   if (nu > 0)
     660             :   {
     661         441 :     GEN u = NULL;
     662             :     while (1)
     663             :     {
     664         444 :       u = lll_block(m, ns+nc, n, ns+nc+ncm, n-1-ncm); if (u) break;
     665           3 :       mprec = precdbl(mprec);
     666             :       /* recompute m0 * u0 to higher prec */
     667           3 :       m = gcharmatnewprec_shallow(gc, mprec);
     668             :     }
     669         441 :     u0 = ZM_mul(u0, u); m = RgM_ZM_mul(m, u);
     670         441 :     if (DEBUGLEVEL>1) err_printf("after LLL reduction (trans block) -> %Ps\n", m);
     671             :   }
     672         532 :   gchar_set_u0(gc, u0); return m;
     673             : }
     674             : 
     675             : /* convert to snf basis of torsion + Z^(r1+2*r2-1) */
     676             : static void
     677         532 : gchar_snfbasis_shallow(GEN gc, GEN rel)
     678             : {
     679             :   long n, r1, r2;
     680             :   GEN nf, cyc, U, Ui;
     681             : 
     682         532 :   nf = gchar_get_nf(gc);
     683         532 :   nf_get_sign(nf, &r1, &r2);
     684         532 :   n = r1+2*r2;
     685             : 
     686         532 :   rel = ZM_hnf(rel);
     687         532 :   if (DEBUGLEVEL>1) err_printf("relations after hnf: %Ps\n", rel);
     688         532 :   cyc = ZM_snf_group(rel, &U, &Ui);
     689         532 :   if (lg(cyc)==1)
     690             :   {
     691         231 :     cyc = zerovec(n-1);
     692         231 :     U = shallowconcat(zeromat(n-1,lg(rel)-1),matid(n-1));
     693         231 :     Ui = shallowtrans(U);
     694             :   }
     695         301 :   else if (n!=1)
     696             :   {
     697         287 :     cyc = shallowconcat(cyc, zerovec(n-1));
     698         287 :     U = shallowmatconcat(mkmat22(U,gen_0,gen_0,matid(n-1)));
     699         287 :     Ui = shallowmatconcat(mkmat22(Ui,gen_0,gen_0,matid(n-1)));
     700             :   }
     701             : 
     702         532 :   rel = shallowtrans(rel);
     703         532 :   U = shallowtrans(U);
     704         532 :   Ui = shallowtrans(Ui);
     705             : 
     706         532 :   gchar_set_nfree(gc, n-1);
     707         532 :   gchar_set_ntors(gc, (lg(cyc)-1) - (n-1));
     708         532 :   gchar_set_cyc(gc, shallowconcat(cyc, real_0(gchar_get_prec(gc))));
     709         532 :   gchar_set_HUUi(gc, rel, U, Ui);
     710         532 : }
     711             : 
     712             : static long
     713        2602 : mextraprec(GEN m_inv, GEN m, GEN gc)
     714             : {
     715        2602 :   return nbits2extraprec(2*maxss(gexpo(m_inv),1) + expu(lg(m))
     716        2602 :           + gexpo(gchar_get_u0(gc)) + 10);
     717             : }
     718             : 
     719             : /* c) transpose inverse + clean rationals.
     720             :    prec = target prec for m^-1,
     721             :    mprec = prec of m */
     722             : static void
     723        1412 : gcharmat_tinverse(GEN gc, GEN m, long prec)
     724             : {
     725             :   GEN m_inv;
     726        1412 :   long k, n, r1, r2, ns, nc, ncm, nm, bitprec, mprec, targetmprec = 0;
     727        1412 :   bitprec = prec2nbits(prec);
     728             : 
     729        1412 :   nf_get_sign(gchar_get_nf(gc), &r1, &r2);
     730        1412 :   n = r1+2*r2;
     731        1412 :   ns = gchar_get_ns(gc);
     732        1412 :   nc = gchar_get_nc(gc);
     733        1412 :   ncm = gchar_get_nalg(gc);
     734        1412 :   nm = ns+nc+n; /* ns + nc + r1 + r2 + r2 */
     735        1412 :   mprec = gprecision(m);
     736             : 
     737             :   while (1)
     738         917 :   {
     739             :     GEN v0, mm;
     740             :     /* insert at column ns+nc+r1+r2, or last column if cm */
     741        2329 :     v0 = vec_v0(nm, ns+nc+1, r1, r2);
     742        2329 :     mm = shallowmatinsert(m, v0, ncm? nm: nm-r2);
     743        2329 :     if (DEBUGLEVEL>1) err_printf("add v0 -> %Ps\n", mm);
     744        2329 :     mm = shallowtrans(mm);
     745        2329 :     m_inv = RgM_inv(mm); /* invert matrix, may need to increase prec */
     746        2329 :     if (m_inv)
     747             :     {
     748        2316 :       if (DEBUGLEVEL>1) err_printf("inverse: %Ps\n",m_inv);
     749        2316 :       m_inv = vecsplice(m_inv, ncm? nm: nm-r2); /* remove v0 */
     750        2316 :       if (DEBUGLEVEL>1) err_printf("v0 removed: %Ps\n", m_inv);
     751        2316 :       m_inv = shallowtrans(m_inv);
     752             :       /* enough precision? */
     753             :       /* |B - A^(-1)| << |B|.|Id-B*A| */
     754        2316 :       if (gexpo(m_inv) + gexpo(gsub(RgM_mul(m_inv, m), gen_1)) + expu(lg(m))
     755        2316 :           <= -bitprec)
     756             :       {
     757             :         /* |A^(-1) - (A+H)^(-1)| << |H|.|A^(-1)|^2 */
     758        1722 :         targetmprec = prec + mextraprec(m_inv,m,gc);
     759        1722 :         if (mprec >= targetmprec) break;
     760             :       }
     761         594 :       else targetmprec = 0;
     762             :     }
     763         917 :     mprec = maxss(precdbl(mprec), targetmprec);
     764         917 :     if (mprec == 2) mprec = 3;
     765         917 :     m = gcharmatnewprec_shallow(gc, mprec); /* m0 * u0 to higher prec */
     766             :   }
     767             :   /* clean rationals */
     768        1412 :   if (nc)
     769             :   { /* reduce mod exponent of the group, TODO reduce on each component */
     770         901 :     GEN zmcyc = locs_get_cyc(gchar_get_zm(gc));
     771         901 :     GEN e = ZV_lcm(zmcyc);
     772             :     long k;
     773        3483 :     for (k = 1; k <= nc; k++)
     774        2582 :       shallow_clean_rat(gel(m_inv, ns+k), 1, nm - 1, /*zmcyc[k]*/e, prec);
     775             :   }
     776        1412 :   if (r2)
     777             :   {
     778        2498 :     for (k = 1; k <= r2; k++)
     779        1544 :       shallow_clean_rat(gel(m_inv,nm-k+1), 1, nm-1, NULL, prec);
     780             :   }
     781        1412 :   if (ncm)
     782             :   {
     783             :     long i, j, e;
     784        1413 :     for (i = 1; i<=r1+r2; i++)
     785        2742 :       for (j = 1; j <= ncm; j++)
     786             :       {
     787        1835 :         e = gexpo(gcoeff(m_inv, ns+nc+j, ns+nc+i));
     788        1835 :         if (e > -20) /* TODO is this bound too permissive? */
     789             :           pari_err_BUG("gcharinit (nonzero entry)"); /*LCOV_EXCL_LINE*/
     790        1835 :         gcoeff(m_inv, ns+nc+j, ns+nc+i) = gen_0;
     791             :       }
     792             :   }
     793        1412 :   if (DEBUGLEVEL>1) err_printf("cyc/cm cleaned: %Ps", shallowtrans(m_inv));
     794             :   /* normalize characters, parameters mod Z */
     795        4799 :   for (k = 1; k <= ns+nc; k++) gel(m_inv, k) = gfrac(gel(m_inv, k));
     796             :   /* increase relative prec of real values */
     797        1412 :   gchar_set_basis(gc, gprec_w(m_inv, prec));
     798        1412 : }
     799             : 
     800             : /* make sure same argument determination is chosen */
     801             : static void
     802        4479 : same_arg(GEN v1, GEN v2, long s1, long s2)
     803             : {
     804        4479 :   long i1, i2, l = lg(v1);
     805       22512 :   for (i1 = s1, i2 = s2; i1 < l; i1++,i2++)
     806             :   {
     807       18033 :     GEN d = grndtoi(gsub(gel(v2,i2),gel(v1,i1)), NULL);
     808       18033 :     if (signe(d)) gel(v1,i1) = gadd(gel(v1,i1), d);
     809             :   }
     810        4479 : }
     811             : 
     812             : static void
     813        4479 : vaffect_shallow(GEN x, long i0, GEN y)
     814             : {
     815             :   long i;
     816       39261 :   for (i = 1; i < lg(y); i++)
     817       34782 :     gel(x, i+i0) = gel(y, i);
     818        4479 : }
     819             : 
     820             : /* recompute matrix with precision increased */
     821             : /* u0 the base change, returns m0 * u0 */
     822             : /* mprec: requested precision for m0 */
     823             : static GEN
     824        1821 : gcharmatnewprec_shallow(GEN gc, long mprec)
     825             : {
     826             :   GEN nf, m0, u0, sfu;
     827             :   long k, l, ns, nc, r1, r2, nfprec, embprec;
     828             : 
     829        1821 :   ns = gchar_get_ns(gc);
     830        1821 :   nc = gchar_get_nc(gc);
     831        1821 :   nf = gchar_get_nf(gc);
     832        1821 :   sfu = gchar_get_sfu(gc);
     833        1821 :   nf_get_sign(nf, &r1, &r2);
     834        1821 :   nfprec = nf_get_prec(gchar_get_nf(gc));
     835             : 
     836        1821 :   m0 = gchar_get_m0(gc);
     837        1821 :   u0 = gchar_get_u0(gc);
     838             : 
     839        1821 :   if (DEBUGLEVEL>3) err_printf("gcharmatnewprec_shallow mprec=%d nfprec=%d\n", mprec, nfprec);
     840             : 
     841        1821 :   embprec = mprec; l = lg(sfu);
     842             :   while(1)
     843             :   { /* recompute the nfembedlogs of s-units and fundamental units */
     844        6300 :     for (k = 1; k < l; k++) /* ns s-units, then r1+r2-1 fundamental units */
     845             :     {
     846        4479 :       GEN emb = nfembedlog(&nf, gel(sfu,k), embprec);
     847        4479 :       if (!emb) break;
     848        4479 :       same_arg(emb, gel(m0,k),  r1+r2, ns+nc+r1+r2);
     849        4479 :       vaffect_shallow(gel(m0, k), ns+nc, emb);
     850             :     }
     851        1821 :     if (k == l) break;
     852           0 :     if (DEBUGLEVEL)
     853           0 :       err_printf("gcharmatnewprec_shallow: increasing embprec %d -> %d\n",
     854             :                  embprec, precdbl(embprec));
     855           0 :     embprec = precdbl(embprec);
     856             :   }
     857        1821 :   gchar_set_nf(gc, nf);
     858        1821 :   gchar_set_nfprec(gc, nfprec);
     859        1821 :   return gmul(m0, u0);
     860             : }
     861             : 
     862             : static void _check_gchar_group(GEN gc, long flag);
     863             : static void
     864       10857 : check_gchar_group(GEN gc) { _check_gchar_group(gc, 0); }
     865             : 
     866             : /* increase prec if needed. FIXME: hardcodes gc[8] and gc[11] */
     867             : /* newprec: requested precision for m_inv */
     868             : static GEN
     869        1645 : gcharnewprec_i(GEN gc, long newprec)
     870             : {
     871             :   long prec, prec0, nfprec, nfprec0, mprec;
     872        1645 :   GEN gc2 = shallowcopy(gc);
     873             : 
     874        1645 :   _check_gchar_group(gc2, 1); /* ignore illegal prec */
     875        1624 :   prec = gchar_get_prec(gc2);
     876        1624 :   nfprec = gchar_get_nfprec(gc2);
     877             : 
     878        1624 :   if (newprec > prec)
     879             :   { /* increase precision */
     880        1118 :     if (DEBUGLEVEL) pari_warn(warnprec,"gcharnewprec",newprec);
     881        1118 :     nfprec += newprec - prec;
     882        1118 :     prec = newprec;
     883        1118 :     gel(gc2, 8) = shallowcopy(gel(gc,8));
     884        1118 :     gmael(gc2, 8, 1) = shallowcopy(gmael(gc, 8, 1));
     885        1118 :     gchar_set_prec(gc2, prec);
     886        1118 :     gchar_set_nfprec(gc2, nfprec);
     887             :   }
     888             : 
     889        1624 :   nfprec0 = nf_get_prec(gchar_get_nf(gc2));
     890        1624 :   if (nfprec0 && nfprec > nfprec0)
     891             :   {
     892         515 :     if (DEBUGLEVEL) pari_warn(warnprec,"gcharnewprec (nf)", nfprec);
     893         515 :     gchar_set_nf(gc2, nfnewprec_shallow(gchar_get_nf(gc2), nfprec));
     894             :   }
     895             : 
     896        1624 :   prec0 = gprecision(gchar_get_basis(gc2));
     897        1624 :   if (prec0 && prec > prec0)
     898             :   {
     899             :     GEN m, cyc;
     900         880 :     if (DEBUGLEVEL) pari_warn(warnprec,"gcharnewprec (minv)", prec);
     901         880 :     gel(gc2, 11) = shallowcopy(gel(gc2, 11));
     902         880 :     mprec = prec + mextraprec(gchar_get_basis(gc), gchar_get_m0(gc), gc);
     903         880 :     m = gcharmatnewprec_shallow(gc2, mprec);
     904         880 :     if (DEBUGLEVEL>2) err_printf("m0*u0 recomputed -> %Ps\n", m);
     905         880 :     gcharmat_tinverse(gc2, m, prec);
     906         880 :     cyc = shallowcopy(gchar_get_cyc(gc2));
     907         880 :     gel(cyc, lg(cyc)-1) = real_0(prec);
     908         880 :     gchar_set_cyc(gc2, cyc);
     909             :   }
     910        1624 :   return gc2;
     911             : }
     912             : 
     913             : /* newprec: requested evalprec */
     914             : GEN
     915        1645 : gcharnewprec(GEN gc, long newprec)
     916             : {
     917        1645 :   pari_sp av = avma;
     918        1645 :   GEN gc2 = gcharnewprec_i(gc, newprec + EXTRAPRECWORD);
     919        1624 :   gchar_set_evalprec(gc2, newprec);
     920        1624 :   return gerepilecopy(av, gc2);
     921             : }
     922             : 
     923             : static void
     924       12488 : check_localstar(GEN x)
     925             : {
     926       12488 :   if (typ(x) != t_VEC || lg(x) != LOCS_LENGTH + 1)
     927           7 :     pari_err_TYPE("char group", x);
     928             :   /* FIXME: check components once settled */
     929       12481 : }
     930             : 
     931             : int
     932        2072 : is_gchar_group(GEN gc)
     933             : {
     934        2072 :   return (typ(gc) == t_VEC
     935        2072 :       &&  lg(gc) == GC_LENGTH + 1
     936         287 :       &&  typ(gel(gc, 8)) == t_VEC
     937         287 :       &&  lg(gel(gc, 8)) == 3
     938         287 :       &&  typ(gmael(gc, 8, 1))  == t_VECSMALL
     939         287 :       &&  typ(gmael(gc, 8, 2))  == t_VECSMALL
     940         287 :       &&  (checkbnf_i(gchar_get_bnf(gc)) != NULL)
     941        4144 :       &&  (checknf_i(gchar_get_nf(gc)) != NULL));
     942             : }
     943             : 
     944             : /* validates the structure format.
     945             :  * Raise error if inconsistent precision, unless flag=1. */
     946             : static void
     947       12502 : _check_gchar_group(GEN gc, long flag)
     948             : {
     949             :   GEN b, bnf, nf;
     950       12502 :   if (typ(gc) != t_VEC || lg(gc) != GC_LENGTH + 1)
     951          14 :     pari_err_TYPE("char group", gc);
     952       12488 :   check_localstar(gchar_get_zm(gc));
     953       12481 :   if (typ(gchar_get_loccyc(gc)) != t_VEC)
     954           0 :     pari_err_TYPE("gchar group (loccyc)", gc);
     955       12481 :   b = gchar_get_basis(gc);
     956       12481 :   if (typ(b) != t_MAT) pari_err_TYPE("gchar group (basis)", gc);
     957       12474 :   bnf = gchar_get_bnf(gc); checkbnf(bnf);
     958       12474 :   nf = gchar_get_nf(gc); checknf(nf);
     959       12474 :   if (!gequal(nf_get_pol(nf),nf_get_pol(bnf_get_nf(bnf))))
     960           0 :     pari_err_TYPE("gchar group (nf != bnf.nf)", gc);
     961       12474 :   if (typ(gel(gc,8)) != t_VEC ||typ(gmael(gc,8,1)) != t_VECSMALL)
     962           7 :     pari_err_TYPE("gchar group (gc[8])", gc);
     963       12467 :   if (!flag)
     964             :   {
     965       10843 :     long prec0 = gprecision(b), nfprec0 = nf_get_prec(nf);
     966       10843 :     if ((prec0 && gchar_get_prec(gc) > prec0)
     967       10836 :         || (nfprec0 && gchar_get_nfprec(gc) > nfprec0))
     968           7 :       pari_err_PREC("gchar group, please call gcharnewprec");
     969             :   }
     970       12460 : }
     971             : 
     972             : /* basis of algebraic characters + cyc vector */
     973             : static GEN
     974          56 : gchar_algebraic_basis(GEN gc)
     975             : {
     976             :   long ntors, nfree, nc, nm, r2, nalg, n0, k;
     977             :   GEN basis, args, m, w, normchar, alg_basis, tors_basis;
     978             : 
     979             :   /* in snf basis */
     980          56 :   ntors = gchar_get_ntors(gc); /* number of torsion generators */
     981          56 :   nfree = gchar_get_nfree(gc);
     982          56 :   nc = ntors + nfree;
     983             :   /* in internal basis */
     984          56 :   n0 = gchar_get_ns(gc) + gchar_get_nc(gc); /* last index of torsion char */
     985          56 :   r2 = gchar_get_r2(gc);
     986          56 :   nm = gchar_get_nm(gc);
     987             :   /* in both */
     988          56 :   nalg = gchar_get_nalg(gc); /* number of generators of free algebraic chars */
     989             : 
     990             :   /* finite order characters have weight 0 */
     991          56 :   tors_basis = matid(ntors);
     992             : 
     993             :   /* the norm is an algebraic character */
     994          56 :   normchar = gneg(col_ei(nc+1,nc+1));
     995             : 
     996             :   /* add sublattice of algebraic */
     997             : 
     998          56 :   if (!nalg)
     999             :   {
    1000          14 :     if (DEBUGLEVEL>2) err_printf("nalg=0\n");
    1001          14 :     return shallowmatconcat(mkvec2(tors_basis,normchar));
    1002             :   }
    1003             : 
    1004             :   /* block of k_s parameters of free algebraic */
    1005          42 :   args = matslice(gchar_get_basis(gc),n0+1,n0+nalg,nm-r2+1,nm);
    1006             : 
    1007          42 :   if (r2 == 1)
    1008             :   {
    1009             :     /* no parity condition */
    1010          14 :     if (DEBUGLEVEL>2) err_printf("r2 = 1 -> args = %Ps\n", args);
    1011          14 :     alg_basis = matid(nalg);
    1012          14 :     w = gel(args,1);
    1013             :   }
    1014             :   else
    1015             :   {
    1016             :     /* parity condition: w + k_s = 0 mod 2 for all s,
    1017             :        ie solve x.K constant mod 2, ie solve
    1018             :        x.K' = 0 mod 2 where K' = [ C-C0 ] (substract first column)
    1019             :      */
    1020             :     /* select block k_s in char parameters and */
    1021          28 :     if (DEBUGLEVEL>2) err_printf("block ks -> %Ps\n", args);
    1022          28 :     m = cgetg(r2, t_MAT);
    1023          77 :     for (k = 1; k < r2; k++)
    1024          49 :       gel(m,k) = gsub(gel(args,k+1),gel(args,1));
    1025          28 :     if (DEBUGLEVEL>2) err_printf("block ks' -> %Ps", m);
    1026          28 :     alg_basis = shallowtrans(gel(matsolvemod(shallowtrans(m),gen_2,gen_0,1),2));
    1027          28 :     if (DEBUGLEVEL>2) err_printf("alg_basis -> %Ps\n", alg_basis);
    1028          28 :     w = gmul(alg_basis, gel(args,1));
    1029          28 :     if (DEBUGLEVEL>2) err_printf("w -> %Ps\n", w);
    1030             :   }
    1031             :   /* add weight to infinite order characters, at position nc+1 */
    1032          42 :   w = gneg(gdivgs(gmodgs(w, 2), 2));
    1033          42 :   alg_basis = shallowmatconcat(mkvec3(alg_basis, zerovec(nfree-nalg), w));
    1034             : 
    1035          42 :   if (lg(tors_basis)>1)
    1036          28 :     basis = shallowmatconcat(mkmat22(tors_basis, gen_0, gen_0, alg_basis));
    1037             :   else
    1038          14 :     basis = alg_basis;
    1039          42 :   return shallowmatconcat(mkvec2(shallowtrans(basis),normchar));
    1040             : }
    1041             : static GEN
    1042         119 : gchar_algebraicnormtype(GEN gc, GEN type)
    1043             : {
    1044         119 :   GEN w = NULL, w1, v;
    1045             :   long i, r1, r2, n;
    1046         119 :   r1 = gchar_get_r1(gc);
    1047         119 :   r2 = gchar_get_r2(gc);
    1048         161 :   for (i=1; i<=r1; i++)
    1049             :   {
    1050          42 :     w1 = addii(gmael(type,i,1),gmael(type,i,2));
    1051          42 :     if (!w) w = w1;
    1052          14 :     else if (!equalii(w,w1)) return NULL;
    1053             :   }
    1054         168 :   for (i=r1+1; i<=r1+r2; i++)
    1055             :   {
    1056         133 :     w1 = gmael(type,i,1);
    1057         133 :     if (!w) w = w1;
    1058          42 :     else if (!equalii(w,w1)) return NULL;
    1059         119 :     if (!equalii(w,gmael(type,i,2))) return NULL;
    1060             :   }
    1061          35 :   n = lg(gchar_get_cyc(gc))-1;
    1062          35 :   v = zerocol(n);
    1063          35 :   gel(v,n) = negi(w);
    1064          35 :   return mkvec(v);
    1065             : }
    1066             : static GEN
    1067         175 : gchar_algebraicoftype(GEN gc, GEN type)
    1068             : {
    1069             :   long i, r1, r2, nalg, n0, nm;
    1070             :   GEN p, q, w, k, matk, chi;
    1071             :   /* in internal basis */
    1072         175 :   n0 = gchar_get_ns(gc) + gchar_get_nc(gc); /* last index of torsion chars */
    1073         175 :   r1 = gchar_get_r1(gc);
    1074         175 :   r2 = gchar_get_r2(gc);
    1075         175 :   nm = gchar_get_nm(gc);
    1076             :   /* in both */
    1077         175 :   nalg = gchar_get_nalg(gc); /* number of generators of free algebraic chars */
    1078             : 
    1079         175 :   if (typ(type)!=t_VEC || lg(type) != r1+r2+1)
    1080          28 :     pari_err_TYPE("gcharalgebraic (1)", type);
    1081         364 :   for (i = 1; i<=r1+r2; i++)
    1082         245 :     if (typ(gel(type,i)) != t_VEC
    1083         238 :         ||lg(gel(type,i)) != 3
    1084         231 :         ||typ(gmael(type,i,1)) != t_INT
    1085         224 :         ||typ(gmael(type,i,2)) != t_INT)
    1086          28 :       pari_err_TYPE("gcharalgebraic (2)", type);
    1087             : 
    1088         119 :   chi = gchar_algebraicnormtype(gc, type);
    1089         119 :   if (chi) return chi;
    1090             : 
    1091          84 :   if (!nalg) return NULL;
    1092          70 :   k = cgetg(r2+1,t_VEC);
    1093          70 :   p = gmael(type, 1, 1);
    1094          70 :   q = gmael(type, 1, 2); w = addii(p, q);
    1095          70 :   gel(k, 1) = subii(q, p);
    1096         105 :   for (i = 2; i <= r2; i++)
    1097             :   {
    1098          49 :     p = gmael(type, i, 1);
    1099          49 :     q = gmael(type, i, 2);
    1100          49 :     if (!equalii(w, addii(p, q))) return NULL;
    1101          35 :     gel(k, i) = subii(q, p);
    1102             :   }
    1103             :   /* block of k_s parameters of free algebraic */
    1104          56 :   matk = matslice(gchar_get_basis(gc),n0+1,n0+nalg,nm-r2+1,nm);
    1105          56 :   chi = inverseimage(shallowtrans(matk),shallowtrans(k));
    1106          56 :   if (lg(chi) == 1) return NULL;
    1107         126 :   for (i=1; i<lg(chi); i++) if (typ(gel(chi,i)) != t_INT) return NULL;
    1108          84 :   chi = mkvec4(zerocol(gchar_get_ntors(gc)), chi,
    1109          42 :                zerocol(gchar_get_nfree(gc) - nalg), gmul2n(negi(w),-1));
    1110          42 :   return mkvec(shallowconcat1(chi));
    1111             : }
    1112             : 
    1113             : GEN
    1114         231 : gcharalgebraic(GEN gc, GEN type)
    1115             : {
    1116         231 :   pari_sp av = avma;
    1117             :   GEN b;
    1118         231 :   check_gchar_group(gc);
    1119         231 :   b = type? gchar_algebraicoftype(gc, type): gchar_algebraic_basis(gc);
    1120         175 :   if (!b) { set_avma(av); return cgetg(1, t_VEC); }
    1121         133 :   return gerepilecopy(av, b);
    1122             : }
    1123             : 
    1124             : /*********************************************************************/
    1125             : /*                                                                   */
    1126             : /*                          CHARACTERS                               */
    1127             : /*                                                                   */
    1128             : /*********************************************************************/
    1129             : /* return chi without (optional) norm component; set *s = to the latter  */
    1130             : static GEN
    1131        8477 : check_gchar_i(GEN chi, long l, GEN *s)
    1132             : {
    1133        8477 :   long i, n = lg(chi)-1;
    1134        8477 :   if (n == l)
    1135             :   { /* norm component */
    1136        2009 :     if (s)
    1137             :     {
    1138        1995 :       *s = gel(chi,l);
    1139        1995 :       switch(typ(*s))
    1140             :       {
    1141        1988 :         case t_INT:
    1142             :         case t_FRAC:
    1143             :         case t_REAL:
    1144        1988 :         case t_COMPLEX: break;
    1145           7 :         default: pari_err_TYPE("check_gchar_i [s]", *s);
    1146             :       }
    1147          14 :     }
    1148        2002 :     chi = vecslice(chi, 1, l-1);
    1149             :   }
    1150             :   else
    1151             :   { /* missing norm component */
    1152        6468 :     if (n != l-1) pari_err_DIM("check_gchar_i [chi]");
    1153        6461 :     if (s) *s = gen_0;
    1154             :   }
    1155       41335 :   for (i = 1; i < l; i++) if (typ(gel(chi,i))!=t_INT)
    1156           7 :     pari_err_TYPE("check_gchar_i [coefficient]", gel(chi,i));
    1157        8456 :   return chi;
    1158             : }
    1159             : 
    1160             : static GEN
    1161        6650 : check_gchar(GEN gc, GEN chi, GEN *s)
    1162             : {
    1163        6650 :   if (typ(chi)!=t_COL) pari_err_TYPE("check_gchar [chi]", chi);
    1164        6629 :   return check_gchar_i(chi, lg(gchar_get_cyc(gc))-1, s);
    1165             : }
    1166             : 
    1167             : static GEN
    1168        1848 : check_gchari(GEN gc, GEN chi, GEN *s)
    1169             : {
    1170        1848 :   if (typ(chi)!=t_VEC) pari_err_TYPE("check_gchari [chi]", chi);
    1171        1848 :   return check_gchar_i(chi, lg(gel(gchar_get_basis(gc),1)), s);
    1172             : }
    1173             : 
    1174             : /* from coordinates on snf basis, return coordinates on internal basis, set
    1175             :  * s to the norm component */
    1176             : static GEN
    1177        6650 : gchar_internal(GEN gc, GEN chi, GEN *s)
    1178             : {
    1179        6650 :   chi = check_gchar(gc, chi, s);
    1180        6608 :   return ZV_ZM_mul(chi, gchar_get_Ui(gc));
    1181             : }
    1182             : 
    1183             : static GEN
    1184       23170 : RgV_frac_inplace(GEN v, long n)
    1185             : {
    1186             :   long i;
    1187       50540 :   for (i = 1; i <= n; i++) gel(v,i) = gfrac(gel(v,i));
    1188       23170 :   return v;
    1189             : }
    1190             : /* from internal basis form, return the R/Z components */
    1191             : static GEN
    1192        8162 : gchari_duallog(GEN gc, GEN chi)
    1193             : {
    1194        8162 :   chi = RgV_RgM_mul(chi, gchar_get_basis(gc)); /* take exponents mod Z */
    1195        8162 :   return RgV_frac_inplace(chi, gchar_get_ns(gc) + gchar_get_nc(gc));
    1196             : }
    1197             : 
    1198             : /* chip has ncp = #zm[1][i].cyc components */
    1199             : static GEN
    1200        1183 : conductor_expo_pr(GEN gens_fil, GEN chip)
    1201             : {
    1202             :   long i;
    1203        2051 :   for (i = lg(gens_fil) - 1; i > 0; i--)
    1204             :   {
    1205        1743 :     GEN gens = gel(gens_fil, i);
    1206             :     long j;
    1207        2863 :     for (j = 1; j < lg(gens); j++)
    1208             :     {
    1209        1995 :       GEN v = gmul(chip, gel(gens, j));
    1210        1995 :       if (denom_i(v) != gen_1) return stoi(i);
    1211             :     }
    1212             :   }
    1213         308 :   return gen_0;
    1214             : }
    1215             : 
    1216             : /* assume chi in log form */
    1217             : static GEN
    1218        1407 : gcharlog_conductor_f(GEN gc, GEN chi, GEN *faN)
    1219             : {
    1220             :   GEN zm, P, E, Lsprk, ufil;
    1221             :   long i, l, ic;
    1222             : 
    1223        1407 :   if (gchar_get_nc(gc) == 0)
    1224             :   {
    1225         252 :     if (faN) *faN = trivial_fact();
    1226         252 :     return gen_1;
    1227             :   }
    1228        1155 :   zm = gchar_get_zm(gc);
    1229        1155 :   Lsprk = locs_get_Lsprk(zm);
    1230        1155 :   ufil = locs_get_Lgenfil(zm);
    1231        1155 :   P = gel(locs_get_famod(zm), 1);
    1232        1155 :   l = lg(Lsprk); E = cgetg(l, t_VEC);
    1233        2338 :   for (i = 1, ic = gchar_get_ns(gc); i < l ; i++)
    1234             :   {
    1235        1183 :     GEN gens = gel(ufil, i), chip;
    1236        1183 :     long ncp = sprk_get_ncp(gel(Lsprk,i));
    1237             : 
    1238        1183 :     chip = vecslice(chi, ic + 1, ic + ncp);
    1239        1183 :     gel(E, i) = conductor_expo_pr(gens, chip);
    1240        1183 :     ic += ncp;
    1241             :   }
    1242        1155 :   if (faN) *faN = famat_remove_trivial(mkmat2(P,E));
    1243        1155 :   return idealfactorback(gchar_get_nf(gc), P, E, 0); /* red = 0 */
    1244             : }
    1245             : 
    1246             : /* ={sigma} s.t. k_sigma = 1 */
    1247             : static GEN
    1248       16415 : gcharlog_conductor_oo(GEN gc, GEN chi)
    1249             : {
    1250       16415 :   long noo, i, n0 = gchar_get_ns(gc) + gchar_get_nc(gc);
    1251             :   GEN k_real, chi_oo, moo, den;
    1252             : 
    1253       16415 :   moo = locs_get_m_infty(gchar_get_zm(gc));
    1254       16415 :   noo = lg(moo)-1;
    1255       16415 :   k_real = vecslice(chi, n0-noo+1, n0);
    1256       16415 :   chi_oo = zerovec(gchar_get_r1(gc));
    1257       17969 :   for (i=1; i<=noo; i++)
    1258             :   {
    1259        1554 :     den = Q_denom(gel(k_real,i));
    1260        1554 :     if (den && !equali1(den)) gel(chi_oo, moo[i]) = gen_1;
    1261             :   }
    1262       16415 :   return chi_oo;
    1263             : }
    1264             : 
    1265             : static GEN
    1266         210 : gchari_conductor(GEN gc, GEN chi)
    1267             : {
    1268         210 :   chi = gchari_duallog(gc, chi);
    1269         210 :   return mkvec2(gcharlog_conductor_f(gc, chi, NULL),
    1270             :                 gcharlog_conductor_oo(gc, chi));
    1271             : }
    1272             : GEN
    1273         210 : gchar_conductor(GEN gc, GEN chi)
    1274             : {
    1275         210 :   pari_sp av = avma;
    1276         210 :   check_gchar_group(gc);
    1277         210 :   return gerepilecopy(av, gchari_conductor(gc, gchar_internal(gc, chi, NULL)));
    1278             : }
    1279             : 
    1280             : int
    1281         175 : gcharisalgebraic(GEN gc, GEN chi, GEN *pq)
    1282             : {
    1283             :   long i, ntors, nfree, n0, nalg, r1, r2;
    1284             :   GEN w, chii, v;
    1285         175 :   pari_sp av = avma;
    1286             : 
    1287         175 :   check_gchar_group(gc);
    1288             :   /* in snf basis */
    1289         175 :   ntors = gchar_get_ntors(gc);
    1290         175 :   nfree = gchar_get_nfree(gc);
    1291             :   /* in internal basis */
    1292         175 :   r1 = gchar_get_r1(gc);
    1293         175 :   r2 = gchar_get_r2(gc);
    1294         175 :   n0 = gchar_get_nm(gc) - r2; /* last index before k_s */
    1295             :   /* in both */
    1296         175 :   nalg = gchar_get_nalg(gc); /* number of generators of free algebraic chars */
    1297             : 
    1298         175 :   chii = gchar_internal(gc, chi, &w);
    1299             :   /* check component are on algebraic generators */
    1300         385 :   for (i = ntors+nalg+1; i <= ntors+nfree; i++)
    1301         252 :     if (!gequal0(gel(chi,i))) return gc_bool(av, 0);
    1302         133 :   chii = gchari_duallog(gc, chii);
    1303             : 
    1304         133 :   if (r1) /* not totally complex: finite order * integral power of norm */
    1305             :   {
    1306          63 :     if (typ(w) != t_INT) return gc_bool(av, 0);
    1307          49 :     w = negi(w);
    1308          49 :     if (pq)
    1309             :     { /* set the infinity type */
    1310          49 :       v = cgetg(r1+r2+1, t_VEC);
    1311         119 :       for (i = 1; i <= r1; i++) gel(v, i) = mkvec2(w, gen_0);
    1312          77 :       for (  ; i <= r1+r2; i++) gel(v, i) = mkvec2(w, w);
    1313             :     }
    1314             :   }
    1315             :   else /* totally complex */
    1316             :   {
    1317             :     int w2;
    1318          70 :     w = gneg(gmul2n(w, 1));
    1319          70 :     if (typ(w) != t_INT) return gc_bool(av, 0);
    1320          63 :     w2 = mpodd(w);
    1321         189 :     for (i = 1; i <= r2; i++) /* need k_s + w = 0 mod 2 for all s */
    1322         133 :       if (mpodd(gel(chii, n0 + i)) != w2) return gc_bool(av, 0);
    1323          56 :     if (pq)
    1324             :     { /* set the infinity type */
    1325          42 :       v = cgetg(r2+1, t_VEC);
    1326         140 :       for (i = 1; i <= r2; i++)
    1327             :       {
    1328          98 :         GEN p = gmul2n(subii(w, gel(chii, n0+i)), -1);
    1329          98 :         gel(v, i) = mkvec2(p, subii(w, p));
    1330             :       }
    1331             :     }
    1332             :   }
    1333         105 :   if (pq) { *pq = gerepilecopy(av, v); av = avma; }
    1334         105 :   return gc_bool(av, 1);
    1335             : }
    1336             : 
    1337             : GEN
    1338         259 : gcharlocal(GEN gc, GEN chi, GEN v, long prec, GEN* pbid)
    1339             : {
    1340         259 :   pari_sp av = avma;
    1341         259 :   GEN nf = gchar_get_nf(gc), s, chiv, logchi;
    1342             : 
    1343         259 :   check_gchar_group(gc);
    1344         259 :   chi = gchar_internal(gc, chi, &s);
    1345         252 :   logchi = gchari_duallog(gc, chi);
    1346         252 :   if (typ(v) == t_INT) /* v infinite */
    1347             :   {
    1348         182 :     long i, r1, r2, tau = itos(v), n0 = gchar_get_ns(gc) + gchar_get_nc(gc);
    1349             :     GEN phi, k;
    1350         182 :     nf_get_sign(nf, &r1, &r2);
    1351         182 :     if (tau <= 0)
    1352           7 :       pari_err_DOMAIN("gcharlocal [index of an infinite place]", "v", "<=", gen_0, v);
    1353         175 :     if (tau > r1+r2)
    1354           7 :       pari_err_DOMAIN("gcharlocal [index of an infinite place]", "v", ">", stoi(r1+r2), v);
    1355         168 :     phi = gel(logchi, n0 + tau);
    1356         168 :     if (tau <= r1) /* v real */
    1357             :     {
    1358          84 :       GEN moo = gel(gchar_get_mod(gc),2);
    1359          84 :       i = zv_search(moo, tau);
    1360          84 :       if (i==0) k = gen_0;
    1361             :       else
    1362             :       {
    1363          21 :         k = gel(logchi, n0 - lg(moo) + 1 + i); /* 0 or 1/2 */
    1364          21 :         if (!gequal0(k)) k = gen_1;
    1365             :       }
    1366             :     }
    1367             :     else /* v complex */
    1368          84 :       k = gel(logchi, n0 + r2 + tau);
    1369         168 :     if (s) phi = gsub(phi, mulcxI(s));
    1370         168 :     chiv = mkvec2(k, phi);
    1371             :   }
    1372             :   else /* v finite */
    1373             :   {
    1374          70 :     GEN P = gchar_get_modP(gc);
    1375             :     long iv;
    1376          70 :     checkprid(v);
    1377          70 :     iv = gen_search(P, v, (void*)cmp_prime_ideal, cmp_nodata);
    1378          70 :     chiv = gchari_eval(gc, NULL, v, 0, logchi, s, prec);
    1379          70 :     if (iv <= 0) chiv = mkvec(chiv);
    1380             :     else
    1381             :     {
    1382          56 :       GEN cyc, w, Lsprk, bid, chip = NULL;
    1383          56 :       long i, ic, l = lg(P);
    1384          56 :       Lsprk = locs_get_Lsprk(gchar_get_zm(gc));
    1385          63 :       for (i = 1, ic = gchar_get_ns(gc); i < l; i++)
    1386             :       {
    1387          63 :         long ncp = sprk_get_ncp(gel(Lsprk,i));
    1388          63 :         if (i == iv) { chip = vecslice(logchi, ic + 1, ic + ncp); break; }
    1389           7 :         ic += ncp;
    1390             :       }
    1391          56 :       if (!chip) pari_err_BUG("gcharlocal (chip not found)");
    1392             :       /* TODO store bid instead of recomputing? */
    1393          56 :       bid = sprk_to_bid(nf, gel(Lsprk,i), nf_INIT);
    1394          56 :       cyc = bid_get_cyc(bid);
    1395          56 :       w = RgV_RgM_mul(chip, gel(bid_get_U(bid),1));
    1396         154 :       for (i = 1; i < lg(w); i++)
    1397          98 :         gel(w,i) = modii(gmul(gel(w,i), gel(cyc,i)), gel(cyc,i));
    1398          56 :       chiv = vec_append(w, chiv);
    1399          56 :       if (pbid) { *pbid = bid; gerepileall(av, 2, &chiv, pbid); return chiv; }
    1400             :     }
    1401             :   }
    1402         210 :   return gerepilecopy(av, chiv);
    1403             : }
    1404             : 
    1405             : 
    1406             : /*******************************************************************/
    1407             : /*                                                                 */
    1408             : /*                EVALUATION OF CHARACTERS                         */
    1409             : /*                                                                 */
    1410             : /*******************************************************************/
    1411             : GEN
    1412        1785 : gcharduallog(GEN gc, GEN chi)
    1413             : {
    1414        1785 :   pari_sp av = avma;
    1415             :   GEN logchi, s;
    1416        1785 :   check_gchar_group(gc);
    1417        1785 :   logchi = gchari_duallog(gc, gchar_internal(gc, chi, &s));
    1418        1785 :   return gerepilecopy(av, shallowconcat1(mkcol2(logchi,s)));
    1419             : }
    1420             : 
    1421             : static GEN
    1422      187845 : gcharisprincipal(GEN gc, GEN x, GEN *palpha)
    1423             : {
    1424      187845 :   GEN t, v, bnf = gchar_get_bnf(gc), DLdata = gchar_get_DLdata(gc);
    1425      187845 :   t = bnfisprincipal0(bnf, x, nf_GENMAT | nf_FORCE); v = gel(t, 1);
    1426      187845 :   *palpha = famat_reduce(famat_mul(nffactorback(bnf, gel(DLdata,2), v),
    1427      187845 :                                    gel(t, 2)));
    1428      187845 :   return ZM_ZC_mul(gel(DLdata,1), v);
    1429             : }
    1430             : 
    1431             : /* complete log of ideal; if logchi != NULL make sure precision is
    1432             :  * sufficient to evaluate gcharlog_eval_raw to attain 'prec' */
    1433             : static GEN
    1434      187845 : gchar_log(GEN gc, GEN x, GEN logchi, long prec)
    1435             : {
    1436      187845 :   GEN zm, v, alpha, arch_log = NULL, zm_log, nf;
    1437             : 
    1438      187845 :   nf = gchar_get_nf(gc);
    1439      187845 :   zm = gchar_get_zm(gc);
    1440      187845 :   v = gcharisprincipal(gc, x, &alpha); /* alpha a GENMAT */
    1441      187845 :   if (DEBUGLEVEL>2) err_printf("v %Ps\n", v);
    1442      187845 :   zm_log = gchar_logm(nf, zm, alpha);
    1443      187845 :   if (DEBUGLEVEL>2) err_printf("zm_log(alpha) %Ps\n", zm_log);
    1444             : 
    1445      187845 :   if (logchi)
    1446             :   { /* check if precision is sufficient, take care of gexpo = -infty */
    1447      184107 :     long e, bit = expu(nf_get_degree(nf) + lg(zm_log)-1) + 3;
    1448      184107 :     e = gexpo(logchi); if (e > 0) bit += e;
    1449      184107 :     e = gexpo(gel(alpha,1)); if (e > 0) bit += e;
    1450      184107 :     prec += nbits2extraprec(bit);
    1451             :   }
    1452             :   for(;;)
    1453             :   {
    1454      187845 :     arch_log = nfembedlog(&nf, alpha, prec);
    1455      187845 :     if (arch_log) break;
    1456           0 :     prec = precdbl(prec);
    1457             :   }
    1458      187845 :   if (DEBUGLEVEL>2) err_printf("arch log %Ps\n", arch_log);
    1459      187845 :   return shallowconcat1(mkvec3(v, gneg(zm_log), gneg(arch_log)));
    1460             : }
    1461             : 
    1462             : /* gp version, with norm component */
    1463             : GEN
    1464          56 : gcharlog(GEN gc, GEN x, long prec)
    1465             : {
    1466          56 :   pari_sp av = avma;
    1467             :   GEN norm;
    1468             : 
    1469          56 :   check_gchar_group(gc);
    1470          56 :   norm = idealnorm(gchar_get_nf(gc), x);
    1471          56 :   norm = mkcomplex(gen_0, gdiv(glog(norm,prec), Pi2n(1,prec)));
    1472          56 :   return gerepilecopy(av, vec_append(gchar_log(gc, x, NULL, prec), norm));
    1473             : }
    1474             : 
    1475             : static GEN
    1476      198611 : gcharlog_eval_raw(GEN logchi, GEN logx)
    1477      198611 : { GEN v = RgV_dotproduct(logchi, logx); return gsub(v, ground(v)); }
    1478             : 
    1479             : /* if flag = 1 -> value in C, flag = 0 -> value in R/Z
    1480             :  * chi in chari format without norm component (given in s)
    1481             :  * if logchi is provided, assume it has enough precision
    1482             :  * TODO check that logchi argument is indeed used correctly by callers */
    1483             : static GEN
    1484      184107 : gchari_eval(GEN gc, GEN chi, GEN x, long flag, GEN logchi, GEN s, long prec)
    1485             : {
    1486             :   GEN z, norm, logx;
    1487      184107 :   if (!logchi)
    1488             :   {
    1489        3920 :     long e, precgc = prec, prec0 = gchar_get_prec(gc);
    1490        3920 :     logchi = gchari_duallog(gc, chi);
    1491        3920 :     e = gexpo(logchi); if (e > 0) precgc += nbits2extraprec(e);
    1492        3920 :     if (precgc > prec0)
    1493             :     {
    1494          14 :       gc = gcharnewprec(gc, precgc);
    1495          14 :       logchi = gchari_duallog(gc, chi);
    1496             :     }
    1497             :   }
    1498      184107 :   logx = gchar_log(gc, x, logchi, prec);
    1499      184107 :   norm = gequal0(s)? NULL: idealnorm(gchar_get_nf(gc), x);
    1500      184107 :   z = gcharlog_eval_raw(logchi, logx);
    1501      184107 :   if (flag)
    1502             :   {
    1503      180117 :     z = expIPiC(gmul2n(z, 1), prec);
    1504      180117 :     if (norm) z = gmul(z, gpow(norm, gneg(s), prec));
    1505             :   }
    1506        3990 :   else if (norm)
    1507         112 :     z = gadd(z, mulcxI(gdiv(gmul(s, glog(norm,prec)), Pi2n(1,prec))));
    1508      184107 :   if (DEBUGLEVEL>1) err_printf("char value %Ps\n", z);
    1509      184107 :   return z;
    1510             : }
    1511             : 
    1512             : GEN
    1513        3934 : gchareval(GEN gc, GEN chi, GEN x, long flag)
    1514             : {
    1515             :   GEN s;
    1516             :   long prec;
    1517        3934 :   pari_sp av = avma;
    1518        3934 :   check_gchar_group(gc);
    1519        3934 :   prec = gchar_get_evalprec(gc);
    1520        3934 :   chi = gchar_internal(gc, chi, &s);
    1521        3920 :   return gerepileupto(av, gchari_eval(gc, chi, x, flag, NULL, s, prec));
    1522             : }
    1523             : 
    1524             : /*******************************************************************/
    1525             : /*                                                                 */
    1526             : /*                         IDENTIFICATION                          */
    1527             : /*                                                                 */
    1528             : /*******************************************************************/
    1529             : static GEN
    1530          98 : col_2ei(long n, long i) { GEN e = zerocol(n); gel(e,i) = gen_2; return e; }
    1531             : static GEN
    1532        3920 : gchar_identify_init(GEN gc, GEN Lv, long prec)
    1533             : {
    1534             :   GEN M, cyc, mult, Lpr, Lk1, Lphi1, Lk2, Llog, eps, U, P, nf, moo, vlogchi;
    1535             :   long beps, r1, r2, d, nmiss, n, npr, nk1, nchi, s, i, j, l, dim, n0, ncol;
    1536             : 
    1537        3920 :   check_gchar_group(gc);
    1538        3899 :   n0 = gchar_get_ns(gc) + gchar_get_nc(gc); /* last index of torsion char */
    1539        3899 :   cyc = gchar_get_cyc(gc);
    1540        3899 :   P = gchar_get_modP(gc);
    1541        3899 :   moo = gel(gchar_get_mod(gc), 2);
    1542        3899 :   nf = gchar_get_nf(gc); nf_get_sign(nf, &r1, &r2);
    1543        3899 :   nchi = lg(cyc)-2; /* ignore norm */
    1544        3899 :   d = r1 + 2*r2;
    1545        3899 :   mult = (nchi >= d)? gel(cyc,1): gen_1;
    1546        3899 :   s = (8*prec2nbits(prec))/10; mult = shifti(mult, s);
    1547        3899 :   beps = -(7*s) / 16; eps = real2n(beps, prec);
    1548        3899 :   l = lg(Lv);
    1549        3899 :   if (lg(gen_sort_uniq(Lv, (void*)cmp_universal, cmp_nodata)) != l)
    1550           7 :     pari_err(e_MISC, "components of Lv must be distinct");
    1551             :   /* log of prime ideals */
    1552        3892 :   Llog = cgetg(l, t_VEC);
    1553             :   /* index in Lchiv corresponding to the places */
    1554        3892 :   Lpr = cgetg(l, t_VECSMALL);
    1555        3892 :   Lk1 = cgetg(l, t_VECSMALL);
    1556        3892 :   Lphi1 = zero_zv(r1);
    1557        3892 :   Lk2 = zero_zv(r2); nmiss = d;
    1558       11564 :   for (i = 1, npr = nk1 = 0; i < l; i++)
    1559        7714 :     if (typ(gel(Lv,i)) == t_INT)
    1560             :     {
    1561        4004 :       long v = itos(gel(Lv,i));
    1562        4004 :       if (v <= 0) pari_err_DOMAIN("gcharidentify", "v", "<=", gen_0, stoi(v));
    1563        3997 :       if (v > r1+r2)
    1564           7 :         pari_err_DOMAIN("gcharidentify", "v", ">", stoi(r1+r2), stoi(v));
    1565        3990 :       if (v <= r1)
    1566             :       { /* don't put in k1 if not in conductor (but keep as phi) */
    1567         196 :         if (zv_search(moo, v)) { nk1++; Lk1[nk1] = i; }
    1568         196 :         Lphi1[v] = i; nmiss--;
    1569             :       }
    1570             :       else
    1571             :       {
    1572        3794 :         Lk2[v-r1] = i; nmiss -= 2;
    1573             :       }
    1574             :     }
    1575             :     else
    1576             :     {
    1577        3710 :       GEN pr = gel(Lv,i);
    1578        3710 :       long lP = lg(P);
    1579        3710 :       if (idealtyp(&pr, NULL) != id_PRIME)
    1580           7 :         pari_err_TYPE("gcharidentify [ideal]", pr);
    1581        3920 :       for (j = 1; j < lP; j++)
    1582         238 :         if (pr_equal(pr, gel(P,j)))
    1583           7 :           pari_err_COPRIME("gcharidentify", pr, gel(gchar_get_mod(gc), 1));
    1584        3682 :       npr++;
    1585        3682 :       Lpr[npr] = i;
    1586        3682 :       gel(Llog,npr) = gchar_log(gc, pr, NULL, prec);
    1587             :     }
    1588        3850 :   setlg(Llog, npr+1); setlg(Lpr, npr+1); setlg(Lk1, nk1+1);
    1589             :   /* build matrix M */
    1590        3850 :   n = npr + nk1; dim = n + d; ncol = n + 1 + nchi;
    1591        3850 :   M = cgetg(ncol+1, t_MAT);
    1592        7532 :   for (j = 1; j <= npr; j++) gel(M,j) = col_ei(dim, j);
    1593        3948 :   for (;  j <= n; j++) gel(M,j) = col_2ei(dim, j);
    1594        3850 :   gel(M,j) = zerocol(dim);
    1595       11585 :   for (i = n+1; i <= n+r1+r2; i++) gcoeff(M,i,j) = eps;
    1596        3850 :   vlogchi = RgM_mul(gchar_get_Ui(gc), gchar_get_basis(gc));
    1597       18858 :   for (j=1; j<=nchi; j++)
    1598             :   {
    1599       15008 :     GEN logchi = RgV_frac_inplace(row(vlogchi, j), n0); /* duallog(e_j) */
    1600       15008 :     GEN chi_oo = gcharlog_conductor_oo(gc, logchi), C = cgetg(dim+1, t_COL);
    1601       29512 :     for (i=1; i<=npr; i++) gel(C,i) = gcharlog_eval_raw(logchi, gel(Llog,i));
    1602       15344 :     for (i=1; i<=nk1; i++) gel(C,npr+i) = gel(chi_oo, itos(gel(Lv,Lk1[i])));
    1603       74774 :     for (i=1; i<=d; i++) gel(C,n+i) = gel(logchi, n0 + i);
    1604       15008 :     gel(M,n+1+j) = C;
    1605             :   }
    1606        4235 :   for (i = 1; i <= r1; i++)
    1607         385 :     if (!Lphi1[i])
    1608             :     {
    1609         189 :       long a = n + i;
    1610         896 :       for (j = 1; j <= ncol; j++) gcoeff(M,a,j) = gmul2n(gcoeff(M,a,j), beps);
    1611             :     }
    1612       11200 :   for (i = 1; i <= r2; i++)
    1613        7350 :     if (!Lk2[i])
    1614             :     {
    1615        3584 :       long a = n + r1 + i, b = a + r2;
    1616       25032 :       for (j = 1; j <= ncol; j++)
    1617             :       {
    1618       21448 :         gcoeff(M,a,j) = gmul2n(gcoeff(M,a,j), beps);
    1619       21448 :         gcoeff(M,b,j) = gmul2n(gcoeff(M,b,j), beps);
    1620             :       }
    1621             :     }
    1622             :   /* rescale and truncate M, then LLL-reduce M */
    1623        3850 :   M = gtrunc(RgM_Rg_mul(M, mult));
    1624        3850 :   U = ZM_lll(M, 0.99, LLL_IM);
    1625        3850 :   M = ZM_mul(M, U);
    1626        3850 :   U = rowslice(U, n + 2, n + 1 + nchi);
    1627        3850 :   return mkvecn(9, M, U, Lpr, Lk1, Lphi1, Lk2, mult, Lv,
    1628             :                    mkvecsmall3(prec, nmiss, beps));
    1629             : }
    1630             : 
    1631             : /* TODO return the distance between the character found and the conditions? */
    1632             : static GEN
    1633        3850 : gchar_identify_i(GEN gc, GEN idinit, GEN Lchiv)
    1634             : {
    1635             :   GEN M, U, Lpr, Lk1, Lphi1, Lk2, mult, cyc, y, x, sumphi, Lv, Norm, nf;
    1636             :   long i, l, r1, r2, beps, npr, nk1, n, nmiss, nnorm, prec;
    1637        3850 :   M = gel(idinit,1);
    1638        3850 :   U = gel(idinit,2);
    1639        3850 :   Lpr = gel(idinit,3);
    1640        3850 :   Lk1 = gel(idinit,4);
    1641        3850 :   Lphi1 = gel(idinit,5);
    1642        3850 :   Lk2 = gel(idinit,6);
    1643        3850 :   mult = gel(idinit,7);
    1644        3850 :   Lv = gel(idinit,8);
    1645        3850 :   prec = gel(idinit,9)[1];
    1646        3850 :   nmiss = gel(idinit,9)[2];
    1647        3850 :   beps = gel(idinit,9)[3];
    1648        3850 :   npr = lg(Lpr)-1;
    1649        3850 :   nk1 = lg(Lk1)-1; n = npr + nk1;
    1650        3850 :   cyc = gchar_get_cyc(gc);
    1651        3850 :   nf = gchar_get_nf(gc);
    1652        3850 :   nf_get_sign(nf, &r1, &r2);
    1653        3850 :   nnorm = 0;
    1654        3850 :   Norm = gen_0;
    1655             : 
    1656        3850 :   l = lg(Lchiv); Lchiv = shallowcopy(Lchiv);
    1657        3850 :   if (lg(Lv) != l) pari_err_DIM("gcharidentify [#Lv != #Lchiv]");
    1658       11438 :   for (i = 1; i < l; i++)
    1659             :   {
    1660        7637 :     GEN t = gel(Lchiv,i), u;
    1661        7637 :     if (typ(gel(Lv,i)) != t_INT)
    1662             :     {
    1663        3682 :       if (typ(t) == t_COMPLEX)
    1664             :       {
    1665           7 :         nnorm++; /* 2 Pi Im(theta) / log N(pr) */
    1666           7 :         Norm = gadd(Norm, gdiv(gmul(Pi2n(1,prec), gel(t,2)),
    1667           7 :                                glog(idealnorm(nf,gel(Lv,i)),prec)));
    1668           7 :         gel(Lchiv,i) = t = gel(t,1);
    1669             :       }
    1670        3682 :       if (!is_real_t(typ(t)))
    1671          14 :         pari_err_TYPE("gcharidentify [character value: should be real or complex]", t);
    1672             :     }
    1673             :     else
    1674             :     {
    1675        3955 :       if (typ(t) != t_VEC)
    1676           7 :         pari_err_TYPE("gcharidentify [character at infinite place: should be a t_VEC]", t);
    1677        3948 :       if (lg(t) != 3)
    1678           7 :         pari_err_DIM("gcharidentify [character at infinite place: should have two components]");
    1679        3941 :       if (typ(gel(t,1)) != t_INT)
    1680           7 :         pari_err_TYPE("gcharidentify [k parameter at infinite place: should be a t_INT]", gel(t,1));
    1681        3934 :       u = gel(t,2);
    1682        3934 :       if (typ(u) == t_COMPLEX)
    1683             :       {
    1684          98 :         nnorm++;
    1685          98 :         Norm = gsub(Norm, gel(u,2)); u = gel(u,1);
    1686          98 :         gel(Lchiv, i) = mkvec2(gel(t,1), u);
    1687             :       }
    1688        3934 :       if (!is_real_t(typ(u)))
    1689           7 :         pari_err_TYPE("gcharidentify [phi parameter at infinite place: should be real or complex]", u);
    1690             :     }
    1691             :   }
    1692             : 
    1693             :   /* construct vector */
    1694        3801 :   y = zerocol(n + r1 + 2*r2); sumphi = gen_0;
    1695        7469 :   for (i=1; i<=npr; i++) gel(y,i) = gel(Lchiv, Lpr[i]);
    1696        3899 :   for (i=1; i<=nk1; i++) gel(y,npr+i) = gmael(Lchiv,Lk1[i],1);
    1697        4186 :   for (i=1; i<=r1; i++)
    1698         385 :     if (Lphi1[i])
    1699             :     {
    1700         196 :       gel(y, n+i) = x =  gmael(Lchiv,Lphi1[i],2);
    1701         196 :       sumphi = gadd(sumphi, x);
    1702             :     }
    1703       11053 :   for (i=1; i<=r2; i++)
    1704        7252 :     if (Lk2[i])
    1705             :     {
    1706        3731 :       long a = n + r1 + i;
    1707        3731 :       gel(y, a + r2) = gmael(Lchiv,Lk2[i],1);
    1708        3731 :       gel(y, a) = x =  gmael(Lchiv,Lk2[i],2);
    1709        3731 :       sumphi = gadd(sumphi, gshift(x,1));
    1710             :     }
    1711        3801 :   if (nmiss)
    1712             :   {
    1713        1918 :     sumphi = gmul2n(gdivgs(sumphi, -nmiss), beps);
    1714        2198 :     for (i = 1; i <= r1; i++) if (!Lphi1[i]) gel(y, n + i) = sumphi;
    1715        5488 :     for (i = 1; i <= r2; i++) if (!Lk2[i])   gel(y, n + r1+i) = sumphi;
    1716             :   }
    1717        3801 :   y = gtrunc(RgC_Rg_mul(y, mult));
    1718             : 
    1719             :   /* find approximation */
    1720        3801 :   x = ZM_ZC_mul(U, RgM_Babai(M, y));
    1721       18613 :   for (i = 1; i < lg(cyc)-1; i++) /* ignore norm */
    1722       14812 :     if (signe(gel(cyc,i))) gel(x,i) = modii(gel(x,i), gel(cyc,i));
    1723        3801 :   if (nnorm > 0) x = vec_append(x, gdivgu(Norm, lg(Lv)-1));
    1724        3801 :   return x;
    1725             : }
    1726             : 
    1727             : /* TODO export the init interface */
    1728             : GEN
    1729        3920 : gchar_identify(GEN gc, GEN Lv, GEN Lchiv, long prec)
    1730             : {
    1731        3920 :   pari_sp av = avma;
    1732        3920 :   GEN idinit = gchar_identify_init(gc, Lv, prec);
    1733        3850 :   return gerepilecopy(av, gchar_identify_i(gc,idinit,Lchiv));
    1734             : }
    1735             : 
    1736             : /*******************************************************************/
    1737             : /*                                                                 */
    1738             : /*                          L FUNCTION                             */
    1739             : /*                                                                 */
    1740             : /*******************************************************************/
    1741             : 
    1742             : /* TODO: could merge with vecan_chigen */
    1743             : 
    1744             : /* return x + yz; y != 0; z = 0,1 "often"; x = 0 "often" */
    1745             : static GEN
    1746     2475592 : gaddmul(GEN x, GEN y, GEN z)
    1747             : {
    1748             :   pari_sp av;
    1749     2475592 :   if (typ(z) == t_INT)
    1750             :   {
    1751     2191245 :     if (!signe(z)) return x;
    1752       13475 :     if (equali1(z)) return gadd(x,y);
    1753             :   }
    1754      293041 :   if (isintzero(x)) return gmul(y,z);
    1755      117362 :   av = avma;
    1756      117362 :   return gerepileupto(av, gadd(x, gmul(y,z)));
    1757             : }
    1758             : 
    1759             : GEN
    1760         630 : vecan_gchar(GEN an, long n, long prec)
    1761             : {
    1762             :   forprime_t iter;
    1763         630 :   GEN gc = gel(an,1), chi = gel(an,2), P = gel(an,3), PZ = gel(an,4);
    1764         630 :   GEN BOUND = stoi(n), v = vec_ei(n, 1);
    1765         630 :   GEN gp = cgetipos(3), nf, chilog, s;
    1766             :   ulong p;
    1767             : 
    1768             :   /* prec increase: 1/n*log(N(pmax)) < log(pmax) */
    1769         630 :   if (DEBUGLEVEL > 1)
    1770           0 :     err_printf("vecan_gchar: need extra prec %ld\n", nbits2extraprec(expu(n)));
    1771         630 :   gc = gcharnewprec(gc, prec + nbits2extraprec(expu(n)));
    1772         630 :   chilog = gchari_duallog(gc, check_gchari(gc, chi, &s));
    1773             : 
    1774         630 :   nf = gchar_get_nf(gc);
    1775             :   /* FIXME: when many log of many primes are computed:
    1776             :      - bnfisprincipal may not be improved
    1777             :      - however we can precompute the logs of generators
    1778             :        for principal part
    1779             :      - when galois, should compute one ideal by orbit.
    1780             :      - when real, clear imaginary part
    1781             :    */
    1782             : 
    1783         630 :   u_forprime_init(&iter, 2, n);
    1784      181405 :   while ((p = u_forprime_next(&iter)))
    1785             :   {
    1786             :     GEN L;
    1787             :     long j;
    1788      180775 :     int check = !umodiu(PZ,p);
    1789      180775 :     gp[2] = p;
    1790      180775 :     L = idealprimedec_limit_norm(nf, gp, BOUND);
    1791      361298 :     for (j = 1; j < lg(L); j++)
    1792             :     {
    1793      180523 :       GEN pr = gel(L, j), ch;
    1794             :       pari_sp av;
    1795             :       ulong k, q;
    1796      180523 :       if (check && gen_search(P, pr, (void*)cmp_prime_ideal, cmp_nodata) > 0)
    1797         420 :         continue;
    1798             :       /* TODO: extract code and use precom sprk? */
    1799      180103 :       av = avma;
    1800      180103 :       ch = gchari_eval(gc, chi, pr, 1, chilog, gen_0, prec);
    1801      180103 :       ch = gerepileupto(av, ch);
    1802      180103 :       q = upr_norm(pr);
    1803      180103 :       gel(v, q) = gadd(gel(v, q), ch);
    1804     2655695 :       for (k = 2*q; k <= (ulong)n; k += q)
    1805     2475592 :         gel(v, k) = gaddmul(gel(v, k), ch, gel(v, k/q));
    1806             :     }
    1807             :   }
    1808             :   /* weight, could consider shifting s at eval instead */
    1809         630 :   if (!gequal0(s))
    1810             :   {
    1811         266 :     GEN ns = dirpowers(n, gneg(s), prec);
    1812             :     long j;
    1813      284781 :     for (j = 1; j <= n; j++)
    1814      284515 :       if (gel(v,j) != gen_0) gel(v, j) = gmul(gel(v,j),gel(ns,j));
    1815             :   }
    1816         630 :   return v;
    1817             : }
    1818             : 
    1819             : GEN
    1820          14 : eulerf_gchar(GEN an, GEN p, long prec)
    1821             : {
    1822          14 :   GEN gc = gel(an,1), chi = gel(an,2), P = gel(an,3), PZ = gel(an,4);
    1823             :   GEN f, L, nf, chilog, s;
    1824             :   int check;
    1825             :   long j, l;
    1826             : 
    1827             :   /* prec increase: 1/n*log(N(pmax)) < log(pmax) */
    1828          14 :   if (DEBUGLEVEL > 1)
    1829           0 :     err_printf("vecan_gchar: need extra prec %ld\n", nbits2extraprec(expi(p)));
    1830          14 :   gc = gcharnewprec(gc, prec + nbits2extraprec(expi(p)));
    1831          14 :   chilog = gchari_duallog(gc, check_gchari(gc, chi, &s));
    1832             : 
    1833          14 :   nf = gchar_get_nf(gc);
    1834          14 :   f = pol_1(0);
    1835          14 :   check = dvdii(PZ, p);
    1836          14 :   L = idealprimedec(nf, p); l = lg(L);
    1837          28 :   for (j = 1; j < l; j++)
    1838             :   {
    1839          14 :     GEN pr = gel(L, j), ch;
    1840          14 :     if (check && gen_search(P, pr, (void*)cmp_prime_ideal, cmp_nodata) > 0)
    1841           0 :       continue;
    1842          14 :     ch =  gchari_eval(gc, chi, pr, 1, chilog, s, prec);
    1843          14 :     f = gmul(f, gsub(gen_1, monomial(ch, pr_get_f(pr), 0)));
    1844             :   }
    1845          14 :   return mkrfrac(gen_1,f);
    1846             : }
    1847             : 
    1848             : static GEN
    1849        1197 : cleanup_vga(GEN vga, long prec)
    1850             : {
    1851             :   GEN ind;
    1852             :   long bitprec, i, l;
    1853        1197 :   if (!prec) return vga; /* already exact */
    1854        1197 :   bitprec = bit_accuracy(prec);
    1855        1197 :   vga = shallowcopy(vga); l = lg(vga);
    1856        4578 :   for (i = 1; i < l; i++)
    1857             :   {
    1858        3381 :     GEN z = gel(vga,i);
    1859        3381 :     if (typ(z) != t_COMPLEX) continue;
    1860        2086 :     if (gexpo(gel(z,2)) < -bitprec+20) gel(vga,i) = gel(z,1);
    1861             :   }
    1862        1197 :   ind = indexsort(imag_i(vga));
    1863        3381 :   for (i = 2; i < l; i++)
    1864             :   {
    1865        2184 :     GEN z = gel(vga,ind[i]), t;
    1866        2184 :     if (typ(z) != t_COMPLEX) continue;
    1867        1456 :     t = imag_i(gel(vga, ind[i-1]));
    1868        1456 :     if (gexpo(gsub(gel(z,2), t)) < -bitprec+20)
    1869         364 :       gel(vga, ind[i]) = mkcomplex(gel(z,1), t);
    1870             :    }
    1871        4578 :   for (i = 1; i < l; i++)
    1872             :   {
    1873        3381 :     GEN z = gel(vga,i);
    1874        3381 :     if (typ(z) != t_COMPLEX) continue;
    1875        2086 :     gel(vga, i) = mkcomplex(gel(z,1), bestappr(gel(z,2), int2n(bitprec/2)));
    1876             :   }
    1877        1197 :   return vga;
    1878             : }
    1879             : 
    1880             : /* TODO: move to lfunutils, use lfunzeta and lfunzetak */
    1881             : GEN
    1882        1204 : gchari_lfun(GEN gc, GEN chi, GEN s0)
    1883             : {
    1884             :   GEN nf, chilog, s, cond_f, cond_oo, vga_r, vga_c, chiw;
    1885             :   GEN v_phi, v_arg, sig, k, NN, faN, P;
    1886             :   long ns, nc, nm, r1, r2;
    1887             : 
    1888        1204 :   nf = gchar_get_nf(gc);
    1889        1204 :   ns = gchar_get_ns(gc);
    1890        1204 :   nc = gchar_get_nc(gc);
    1891        1204 :   nm = gchar_get_nm(gc);
    1892        1204 :   nf_get_sign(nf, &r1, &r2);
    1893        1204 :   chi = check_gchari(gc, chi, &s);
    1894        1204 :   chilog = gchari_duallog(gc, chi);
    1895        1204 :   s = gadd(s0,s); chiw = shallowconcat(chi, s);
    1896        1204 :   if (!gequal0(imag_i(s)))
    1897           7 :     pari_err_IMPL("lfun for gchar with imaginary norm component");
    1898        1197 :   cond_f =  gcharlog_conductor_f(gc, chilog, &faN);
    1899        1197 :   P = gel(faN, 1); /* prime ideals dividing cond(chi) */
    1900        1197 :   cond_oo =  gcharlog_conductor_oo(gc, chilog);
    1901             : 
    1902        1197 :   NN = mulii(idealnorm(nf, cond_f), absi_shallow(nf_get_disc(nf)));
    1903             :   /* FIXME: shift by s */
    1904        1197 :   if (equali1(NN)) return lfuncreate(gen_1);
    1905        1197 :   if (ZV_equal0(chi)) return lfuncreate(nf);
    1906             : 
    1907             :   /* vga_r = vector(r1,k,I*c[ns+nc+k]-s + cond_oo[k]);
    1908             :    * vga_c = vector(r2,k,I*c[ns+nc+r1+k]+abs(c[ns+nc+r1+r2+k])/2-s) */
    1909        1197 :   v_phi = gmul(vecslice(chilog, ns+nc+1, ns+nc+r1+r2), gen_I());
    1910        1197 :   v_arg = gdivgs(gabs(vecslice(chilog,ns+nc+r1+r2+1,nm),BITS_IN_LONG), 2);
    1911        1197 :   vga_r = gadd(vecslice(v_phi, 1, r1), cond_oo);
    1912        1197 :   vga_c = gadd(vecslice(v_phi, r1+1, r1+r2), v_arg);
    1913        1197 :   sig = shallowconcat1(mkvec3(vga_r,vga_c,gadd(vga_c,const_vec(r2,gen_1))));
    1914             :   /* TODO: remove cleanup when gammamellinv takes ldata*/
    1915        1197 :   sig = cleanup_vga(sig, gchar_get_prec(gc));
    1916        1197 :   k = gen_1;
    1917        1197 :   if (!gequal0(s))
    1918             :   {
    1919             :     long j;
    1920        1484 :     for (j = 1; j < lg(sig); j++) gel(sig, j) = gadd(gel(sig, j), s);
    1921         413 :     k = gsub(k, gmulgs(s,2));
    1922             :   }
    1923             : 
    1924             : #define tag(x,t)  mkvec2(mkvecsmall(t),x)
    1925        1197 :   return mkvecn(6, tag(mkvec4(gc, chiw, P, prV_lcm_capZ(P)), t_LFUN_HECKE),
    1926             :                 gen_1, sig, k, NN, gen_0);
    1927             : }
    1928             : 
    1929             : GEN
    1930         287 : lfungchar(GEN gc, GEN chi)
    1931             : {
    1932         287 :   pari_sp av = avma;
    1933             :   GEN s;
    1934         287 :   check_gchar_group(gc);
    1935         287 :   chi = gchar_internal(gc, chi, &s);
    1936         273 :   return gerepilecopy(av, gchari_lfun(gc, chi, s));
    1937             : }

Generated by: LCOV version 1.13