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

Generated by: LCOV version 1.14