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 - hgm.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.18.1 lcov report (development 30005-fc14bb602a) Lines: 1216 1255 96.9 %
Date: 2025-02-18 09:22:46 Functions: 139 139 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2018  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : /****************************************************************/
      15             : /*                   HYPERGEOMETRIC MOTIVES                     */
      16             : /****************************************************************/
      17             : #include "pari.h"
      18             : #include "paripriv.h"
      19             : 
      20             : #define DEBUGLEVEL DEBUGLEVEL_hgm
      21             : static THREAD GEN HGMDATA2, HGMDATA3;
      22             : 
      23             : void
      24      337970 : pari_init_hgm(void) { HGMDATA2 = HGMDATA3 = NULL; }
      25             : 
      26             : void
      27      333120 : pari_close_hgm(void)
      28      333120 : { guncloneNULL_deep(HGMDATA2); guncloneNULL_deep(HGMDATA3); }
      29             : 
      30             : /* For now: only HGM defined over Q. Can be given as
      31             :  * [alpha] two lists (alpha_i), (beta_j) of the same length of
      32             :  *      rational numbers in [0,1[ with the cyclotomic property.
      33             :  * [cyclo] two lists (d_i), (e_j) of positive integers with
      34             :  *      sum_i phi(d_i) = sum_j phi(e_j).
      35             :  * [gamma] a list of (g_i) such that P=prod_i(x^i-1)^{g_i}
      36             :  *
      37             :  * Initialization computes all HGM data independent of t (hgm) */
      38             : /***************************************************************/
      39             : /*            PART I: functions independent of t, p            */
      40             : /***************************************************************/
      41             : 
      42             : static GEN
      43        1064 : RgV_addhalf(GEN x) { pari_APPLY_same(gadd(ghalf, gel(x,i))); }
      44             : static GEN
      45      979433 : RgV_mirror(GEN x) { pari_APPLY_same(gsubsg(1, gel(x,i))); }
      46             : 
      47             : static GEN
      48      104146 : hodge(GEN val, GEN vbe, long *pTT)
      49             : {
      50      104146 :   long r = lg(val) - 1, T, s, k;
      51      104146 :   GEN H, w = indexsort( shallowconcat(val, RgV_mirror(vbe)) );
      52     1854720 :   for (k = 1, s = T = 0; k <= 2*r; k++)
      53     1750574 :     if (w[k] > r) { if (--s < T) T = s; } else s++;
      54      104146 :   H = zero_zv(r + 1 - T);
      55     1854720 :   for (k = 1, s = 0; k <= 2*r; k++) /* (x^-T * Hodge) += x^s, s - T >= 0 */
      56     1750574 :     if (w[k] > r) s--; else { H[s + 1 - T]++; s++; }
      57      104146 :   *pTT = -T; return zv_to_zx(H, evalvarn(0));
      58             : }
      59             : 
      60             : /* Conversion from [val] to cyclotomic factorization: t_VECSMALL E encodes
      61             :  * prod_i polcyclo(i)^E[i] */
      62             : static GEN
      63        1414 : al2cyE(GEN v)
      64             : {
      65             :   GEN dv, E, F, w;
      66        1414 :   long N, i, j, l = lg(v);
      67        1414 :   if (l == 1) return cgetg(1, t_VECSMALL);
      68        1246 :   w = Q_remove_denom(v, &dv);
      69        1246 :   if (!dv) return mkvecsmall(l-1);
      70        1001 :   N = itou(dv); w = ZV_to_Flv(w, N); vecsmall_sort(w);
      71        1001 :   E = zero_zv(N); /* to record polcyclo(i)^E[i] */
      72        1001 :   F = cgetg(l, t_VECSMALL); /* to check input */
      73        4949 :   for (i = j = 1; i < l; i++)
      74             :   {
      75        3948 :     long k, m = w[i];
      76        3948 :     if (!m) { E[1]++; F[j++] = 0; }
      77        3563 :     else if (N % m == 0)
      78             :     {
      79        1904 :       long d = N/m, km;
      80        1904 :       E[d]++;
      81        9037 :       for (k = 1, km = m; k <= d; k++, km += m)
      82        7133 :         if (ugcd(k,d) == 1) F[j++] = km;
      83             :     }
      84             :   }
      85        1001 :   setlg(F,j); vecsmall_sort(F);
      86        1001 :   return gequal(w,F)? E: NULL;
      87             : }
      88             : 
      89             : static int
      90         707 : cyE_intersect(GEN a, GEN b)
      91             : {
      92         707 :   long i, l = minss(lg(a), lg(b));
      93        2695 :   for (i = 1; i < l; i++) if (a[i] && b[i]) return 1;
      94         700 :   return 0;
      95             : }
      96             : static GEN
      97         707 : get_CYCLOE(GEN val, GEN vbe)
      98             : {
      99         707 :   GEN Ea = al2cyE(val), Eb = al2cyE(vbe);
     100         707 :   if (!Ea || !Eb || cyE_intersect(Ea, Eb))
     101           7 :     pari_err_TYPE("hgminit [not a Q motive]", mkvec2(val,vbe));
     102         700 :   return mkvec2(Ea, Eb);
     103             : }
     104             : 
     105             : /* R[n/d] += mu(d) * r, for all d | n */
     106             : static void
     107       10549 : moebiusadd(GEN R, long n, long r)
     108             : {
     109       10549 :   if (r) {
     110        1785 :     GEN D = divisorsu_moebius(gel(factoru(n), 1));
     111        1785 :     long j, l = lg(D);
     112        1785 :     R[n] += r; /* d = 1 */
     113        3696 :     for (j = 2; j < l; j++) { long d = D[j]; R[n / labs(d)] += d < 0? -r: r; }
     114             :   }
     115       10549 : }
     116             : /* Use Phi_n(x) = Prod_{d|n} (x^d-1)^mu(n/d); Ea/Eb in terms of (x^i-1) */
     117             : static GEN
     118         700 : get_VPOLGA(GEN E)
     119             : {
     120         700 :   GEN Ea = gel(E,1), Eb = gel(E,2), R;
     121         700 :   long i, la = lg(Ea), lb = lg(Eb), n = maxuu(la, lb) - 1;
     122             :   pari_sp av;
     123         700 :   R = zero_zv(n); av = avma;
     124        7770 :   for (i = 1; i < la; i++) moebiusadd(R, i, Ea[i]);
     125        4179 :   for (i = 1; i < lb; i++) moebiusadd(R, i, -Eb[i]);
     126        5733 :   for (i = n; i; i--) if (R[i]) break;
     127         700 :   setlg(R,i+1); set_avma(av); return R;
     128             : }
     129             : 
     130             : /* disc(polcyclo(n)) modulo squares */
     131             : static GEN
     132         126 : cyclodiscmodsq(ulong n)
     133             : {
     134             :   long e, s;
     135             :   ulong p;
     136             : 
     137         126 :   if (n <= 2) return gen_1;
     138         126 :   if ((n & 3) == 2) n >>= 1;
     139         126 :   e = uisprimepower(n, &p);
     140         126 :   if (!e) return gen_1;
     141         112 :   if (p == 2) return e == 2? gen_m1: gen_1;
     142          63 :   s = odd(e)? p: 1;
     143          63 :   if ((p & 3) == 3) s = -s;
     144          63 :   return stoi(s);
     145             : }
     146             : /* E encodes prod polcyclo(i)^E(i) */
     147             : static GEN
     148         231 : discprod(GEN E)
     149             : {
     150         231 :   long i, j, l = lg(E);
     151         231 :   GEN P = cgetg(l, t_VEC);
     152        1148 :   for (i = 3, j = 1; i < l; i++) /* can skip polcyclo(1 or 2) */
     153         917 :     if (odd(E[i])) gel(P,j++) = cyclodiscmodsq(i);
     154         231 :   setlg(P,j); return ZV_prod(P);
     155             : }
     156             : 
     157             : static GEN
     158        1218 : QV_normalize(GEN a)
     159             : {
     160        1218 :   long l = lg(a), i;
     161        1218 :   GEN v = cgetg(l, t_VEC);
     162        5978 :   for (i = 1; i < l; i++)
     163             :   {
     164        4774 :     GEN c = gel(a, i);
     165        4774 :     if (!is_rational_t(typ(c)))
     166          14 :       pari_err_TYPE("hgminit [not rational params]",c);
     167        4760 :     gel(v, i) = gfrac(c);
     168             :   }
     169        1204 :   return sort(v);
     170             : }
     171             : 
     172             : static GEN
     173        1365 : mangoldtexp(long n)
     174             : {
     175        1365 :   GEN D = divisorsu_moebius(gel(factoru(n),1)), P = gen_1;
     176        1365 :   long i, l = lg(D);
     177        4151 :   for (i = 1; i < l; i++)
     178             :   { /* prod (n/d) ^ (n/d mu(d)) */
     179        2786 :     long mud = D[i], nd = n / mud; /* mud = mu(d) * d */
     180        2786 :     P = gmul(P, powis(utoi(labs(nd)), nd));
     181             :   }
     182        1365 :   return P;
     183             : }
     184             : static GEN
     185         938 : E2exp(GEN E)
     186             : {
     187         938 :   long n, r, l = lg(E);
     188         938 :   GEN P = gen_1;
     189        9877 :   for (n = 1; n < l; n++)
     190        8939 :     if ((r = E[n])) P = gmul(P, gpowgs(mangoldtexp(n), r));
     191         938 :   return P;
     192             : }
     193             : static long
     194         973 : get_b1(GEN CYCLOE) { return gel(CYCLOE,2)[1]; }
     195             : 
     196             : static GEN get_u(GEN al, GEN be, GEN CYCLOE, GEN VPOLGA, long DEG, long WT);
     197             : /* [a,b] = alpha, beta */
     198             : static GEN
     199         511 : initab(GEN a, GEN b)
     200             : {
     201             :   GEN VALNUM, VALDEN, VBENUM, VBEDEN, CYCLOE, SIGNPAR;
     202             :   GEN BAD, HODGE, VPOLGA, res;
     203         511 :   long j, WT, TT, OFFMPOL, l = lg(a), DEG = l-1, SWAP = 0;
     204             : 
     205         511 :   if (lg(b) != l) pari_err_TYPE("hgminit [#al != #be]", mkvec2(a,b));
     206         511 :   if (l == 1) pari_err_TYPE("hgminit [#al = 0]", mkvec2(a,b));
     207         497 :   a = QV_normalize(a);
     208         490 :   b = QV_normalize(b);
     209         483 :   if (isintzero(gel(a,1)))
     210             :   {
     211          49 :     if (isintzero(gel(b,1)))
     212           7 :       pari_err_TYPE("hgminit [nonempty intersection]", mkvec2(a,b));
     213          42 :     swap(a, b); SWAP = 1;
     214             :   }
     215         476 :   VALNUM = cgetg(l, t_VECSMALL);
     216         476 :   VALDEN = cgetg(l, t_VECSMALL);
     217         476 :   VBENUM = cgetg(l, t_VECSMALL);
     218         476 :   VBEDEN = cgetg(l, t_VECSMALL);
     219        2478 :   for (j = 1; j < l; j++)
     220             :   {
     221        2002 :     Qtoss(gel(a,j), &VALNUM[j], &VALDEN[j]);
     222        2002 :     Qtoss(gel(b,j), &VBENUM[j], &VBEDEN[j]);
     223             :   }
     224         476 :   BAD = negi(lcmii(ZV_lcm(zv_to_ZV(VALDEN)), ZV_lcm(zv_to_ZV(VBEDEN))));
     225         476 :   HODGE = hodge(a, b, &TT);
     226         476 :   WT = degpol(HODGE);
     227         476 :   CYCLOE = get_CYCLOE(a, b);
     228         469 :   VPOLGA = get_VPOLGA(CYCLOE);
     229         469 :   SIGNPAR = odd(WT)? gen_1: discprod(gel(CYCLOE, odd(DEG)? 2: 1));
     230         469 :   OFFMPOL = (get_b1(CYCLOE) - zv_sum(VPOLGA)) >> 1;
     231         469 :   res = cgetg(13, t_VEC);
     232         469 :   gel(res, 1) = a;
     233         469 :   gel(res, 2) = b;
     234         469 :   gel(res, 3) = CYCLOE;
     235         469 :   gel(res, 4) = VBENUM;
     236         469 :   gel(res, 5) = VBEDEN;
     237         469 :   gel(res, 6) = gdiv(E2exp(gel(CYCLOE,1)), E2exp(gel(CYCLOE,2))); /*MVALUE*/
     238         469 :   gel(res, 7) = VPOLGA;
     239         469 :   gel(res, 8) = BAD;
     240         469 :   gel(res, 9) = HODGE;
     241         469 :   gel(res, 10) = get_u(a, b, CYCLOE, VPOLGA, DEG, WT);
     242         469 :   gel(res, 11) = SIGNPAR;
     243         469 :   gel(res, 12) = mkvecsmall3(OFFMPOL, TT, SWAP);
     244         469 :   return res;
     245             : }
     246             : static int
     247        1295 : checkhgm(GEN v)
     248             : {
     249        1274 :   return typ(v) == t_VEC && lg(v) == 13
     250        2569 :       && typ(gel(v,12)) == t_VECSMALL && lg(gel(v,12)) == 4;
     251             : }
     252             : 
     253             : /* accessors */
     254      141822 : static GEN hgm_get_VAL(GEN v) { return gel(v, 1); }
     255          35 : static GEN hgm_get_VBE(GEN v) { return gel(v, 2); }
     256       12446 : static GEN hgm_get_CYCLOE(GEN v) { return gel(v, 3); }
     257       37597 : static GEN hgm_get_VBENUM(GEN v) { return gel(v, 4); }
     258       37597 : static GEN hgm_get_VBEDEN(GEN v) { return gel(v, 5); }
     259         693 : static GEN hgm_get_MVALUE(GEN v) { return gel(v, 6); }
     260      698738 : static GEN hgm_get_VPOLGA(GEN v) { return gel(v, 7); }
     261       37844 : static GEN hgm_get_BAD(GEN v) { return gel(v, 8); }
     262      109796 : static GEN hgm_get_HODGE(GEN v) { return gel(v, 9); }
     263        2940 : static GEN hgm_get_U(GEN v) { return gmael(v, 10,1); }
     264          63 : static GEN hgm_get_U0(GEN v) { return gmael(v, 10,2); }
     265       28628 : static GEN hgm_get_SIGNPAR(GEN v) { return gel(v, 11); }
     266      141788 : static long hgm_get_DEG(GEN v) { return lg(hgm_get_VAL(v))-1; }
     267      652497 : static long hgm_get_OFFMPOL(GEN v) { return gel(v, 12)[1]; }
     268       87831 : static long hgm_get_TT(GEN v) { return gel(v, 12)[2]; }
     269         686 : static long hgm_get_SWAP(GEN v) { return gel(v, 12)[3]; }
     270      109432 : static long hgm_get_WT(GEN v) { return degpol(hgm_get_HODGE(v)); }
     271             : 
     272             : /***************************************************************/
     273             : /*                 PART II: p-adic functions                   */
     274             : /***************************************************************/
     275             : 
     276             : /* Compute p^k u_{pk+r} for 0 <= r < p and 0 <= k < ell */
     277             : static GEN
     278        2155 : compu(long ell, long p, long D)
     279             : {
     280        2155 :   GEN v = cgetg(p + 1, t_VEC), w = cgetg(ell + 2, t_VEC);
     281        2155 :   long s, t, k, i, L = ell * p;
     282             : 
     283        2155 :   gel(v,1) = gel(v,2) = cvtop(gen_1, stoi(p), D);
     284       59335 :   for (i = 2; i < p; i++) gel(v, i+1) = gdivgu(gel(v, i), i);
     285        2155 :   gel(w, 1) = shallowcopy(v);
     286      260843 :   for (k = s = 1, t = p; k <= L; k++)
     287             :   {
     288      258700 :     gel(v, s) = gdivgu(gadd(gel(v, t), gel(v, s)), k + p - 1);
     289      258687 :     if (++s > p) { s = 1; gel(w, k/p + 1) = shallowcopy(v); }
     290      258688 :     if (++t > p) t = 1;
     291             :   }
     292        2143 :   return w;
     293             : }
     294             : 
     295             : /* [binom(x + k - 1, k) k! p^k, k = 1..N], x a t_PADIC */
     296             : static GEN
     297      116951 : binomfact(GEN x, long N)
     298             : {
     299      116951 :   GEN z, v = cgetg(N+1, t_VEC);
     300             :   long i;
     301      116953 :   gel(v,1) = z = leafcopy(x);
     302      116952 :   setvalp(z, valp(z)+1);
     303      625655 :   for (i = 1; i < N; i++)
     304             :   {
     305      508708 :     gel(v, i+1) = z = gmul(z, gaddgs(x, i));
     306      508679 :     setvalp(z, valp(z) + 1);
     307             :   }
     308      116947 :   return v;
     309             : }
     310             : /* gamma(m /p^f-1) mod p^dfp, m = m0 (mod p), x = (m\p + m0*p^(f-1)) / p^f-1 */
     311             : static ulong
     312      116951 : sumbinom(GEN v, GEN x, long m0, long N, ulong pD)
     313             : {
     314      116951 :   pari_sp av = avma;
     315      116951 :   GEN C = binomfact(x, N), S = gmael(v, 1, m0+1); /* k = 0 */
     316             :   long k;
     317      742572 :   for (k = 1; k <= N; k++) S = gadd(S, gmul(gel(C,k), gmael(v, k+1, m0+1)));
     318      116951 :   return gc_ulong(av, Rg_to_Fl(S, pD));
     319             : }
     320             : /* x * y mod p^D, y a p-unit */
     321             : static GEN
     322      203899 : umultop(ulong x, ulong y, GEN gp, GEN gpD, long D)
     323             : {
     324      203899 :   ulong pD = gpD[2], xy;
     325             :   long v;
     326      203899 :   if (!x) return zeropadic_shallow(gp, D);
     327      203899 :   v = u_lvalrem(x, gp[2], &x);
     328      203899 :   if (x >= pD) x %= pD;
     329      203899 :   xy = Fl_mul(x, y, pD);
     330      203899 :   retmkpadic(utoi(xy), gp, gpD, v, D);
     331             : }
     332             : 
     333             : /* Compute all gamma(m/(p^f-1)) mod p^D for 0 <= m < p^f-1 */
     334             : static GEN
     335        2162 : precomp(ulong p, long f, long D)
     336             : {
     337             :   pari_sp av, av2;
     338             :   GEN vall, gp, gpD, vres, x, l2;
     339             :   ulong m, m0, m1, M, q, v, vpo, vva, vpowv, pf1;
     340        2162 :   ulong N = ceildivuu(D * p, p-1), pf = upowuu(p,f), pD = upowuu(p, D), pM;
     341        2162 :   ulong ga12, S, iQ, Q = pf - 1;
     342             : 
     343        2162 :   if (!pD || !pf) pari_err_OVERFLOW("hgmprecomp");
     344        2155 :   iQ = D <= f? pD - 1: Fl_inv(Q, pD);
     345        2155 :   vres = cgetg(pf, t_VECSMALL); vres[1] = 1; av = avma;
     346        2155 :   gp = utoipos(p); gpD = utoipos(pD);
     347        2155 :   vall = compu(N, p, D);
     348        2155 :   if (p == 2)
     349             :   {
     350         133 :     if (f == 1) return gc_const(av, vres);
     351          14 :     av2 = avma;
     352          70 :     for (m = 1; m < Q; m++, set_avma(av2))
     353             :     {
     354          56 :       m0 = m&1, m1 = m >> 1;
     355          56 :       x = umultop(m1 + (m0 << (f - 1)), iQ, gp, gpD, D);
     356          56 :       vres[m+1] = sumbinom(vall, x, m0, N, pD);
     357             :     }
     358          14 :     return gc_const(av, vres);
     359             :   }
     360        2022 :   q = Q >> 1; pf1 = pf / p;
     361        2022 :   vva = z_lval(Q, 2); vpowv = 1 << vva;
     362        2022 :   av2 = avma;
     363      104481 :   for (m = 1; m <= q; m += 2, set_avma(av2))
     364             :   {
     365      102459 :     if (f > 1) { m0 = m%p; m1 = m / p; M = m1 + m0 * pf1; } else M = m0 = m;
     366      102459 :     x = umultop(M, iQ, gp, gpD, D);
     367      102458 :     vres[m+1] = S = sumbinom(vall, x, m0, N, pD);
     368      102459 :     if (!odd(m0)) S = pD - S;
     369      102459 :     vres[pf - m] = Fl_inv(S, pD);
     370             :   }
     371       16460 :   for (m = vpowv; m <= q; m += 2*vpowv, set_avma(av2))
     372             :   {
     373       14438 :     if (f > 1) { m0 = m%p; m1 = m / p; M = m1 + m0 * pf1; } else M = m0 = m;
     374       14438 :     x = umultop(M, iQ, gp, gpD, D);
     375       14438 :     vres[m+1] = S = sumbinom(vall, x, m0, N, pD);
     376       14438 :     if (!odd(m0)) S = pD - S;
     377       14438 :     vres[pf - m] = Fl_inv(S, pD);
     378             :   }
     379        2022 :   l2 = gmulgu(Qp_log(cvtop(gen_2, gp, D)), p - 1);
     380        2022 :   ga12 = Fl_inv(Rg_to_Fl(Qp_gamma(cvtop(ghalf, gp, D)), pD), pD);
     381        2022 :   pM = maxuu(pD, pf); av2 = avma;
     382       10637 :   for (v = 1, vpo = 2; vpo < Q; v++, vpo <<= 1)
     383             :   {
     384        8615 :     if (v == vva) continue;
     385       93797 :     for (m = vpo; m <= q; m += 2*vpo, set_avma(av2))
     386             :     {
     387       86947 :       ulong ms2 = (m >> 1) + 1;
     388       86947 :       ulong t = Fl_mul(Fl_mul(vres[ms2], vres[ms2+q], pD), ga12, pD);
     389             :       GEN mp;
     390             :       /* - m  = -m1 * p + m0, 1 <= m0 <= p */
     391       86946 :       m0 = p - m%p; m1 = m/p + 1; M = pM + m1 - m0 * pf1;
     392       86946 :       mp = umultop(M, iQ, gp, gpD, D);
     393       86945 :       S = Rg_to_Fl(gmul(Qp_exp(gmul(mp, l2)), shifti(utoi(t), m0-1)), pD);
     394       86944 :       vres[m+1] = S;
     395       86944 :       if (odd(m0)) S = pD - S;
     396       86944 :       vres[pf - m] = Fl_inv(S, pD);
     397             :     }
     398             :   }
     399        2022 :   return gc_const(av, vres);
     400             : }
     401             : 
     402             : static long
     403        3661 : get_pad(ulong p)
     404             : {
     405        3661 :   switch(p) {
     406         182 :     case 2: return 18;
     407         399 :     case 3: return 11;
     408         322 :     case 5: return  8;
     409             :   }
     410        2758 :   if (p <= 37) return 6;
     411         819 :   if (p <= 251) return 4;
     412           7 :   if (p <= 65521) return 2;
     413           7 :   return 1;
     414             : }
     415             : 
     416             : static GEN
     417         609 : Flv_red(GEN z, ulong p)
     418             : {
     419         609 :   long i, l = lg(z);
     420         609 :   GEN x = cgetg(l, t_VECSMALL);
     421      248780 :   for (i=1; i<l; i++) x[i] = uel(z,i)%p;
     422         609 :   return x;
     423             : }
     424             : 
     425             : static GEN
     426          97 : block0(long l)
     427             : {
     428          97 :   GEN v = cgetg_block(l, t_VEC);
     429             :   long i;
     430       90066 :   for (i = 1; i < l; i++) gel(v,i) = gen_0;
     431          97 :   return v;
     432             : }
     433             : 
     434             : /* f > 0, p prime, need result mod p^dfp */
     435             : static GEN
     436        3668 : doprecomp(ulong p, long f, long dfp)
     437             : {
     438             :   GEN t, T;
     439             :   long F; /* can compute cache mod p^F < 2^BIL */
     440             :   ulong m, n, q;
     441        3668 :   if (f > 3 || (F = get_pad(p)) < dfp) return precomp(p, f, dfp);
     442        3633 :   if (f == 3)
     443             :   {
     444          35 :     T = HGMDATA3; if (!T) T = HGMDATA3 = block0(100);
     445          35 :     if (p >= (ulong)lg(T)) return precomp(p, f, dfp);
     446          35 :     t = gel(T,p);
     447          35 :     if (typ(t) == t_INT) t = gel(T,p) = gclone(precomp(p, f, F));
     448          35 :     return (F == dfp)? t: Flv_red(t, upowuu(p, dfp));
     449             :   }
     450             :   /* f <= 2, dfp <= F */
     451        3598 :   T = HGMDATA2; if (!T) T = HGMDATA2 = block0(1000);
     452        3598 :   if (p >= (ulong)lg(T)) return precomp(p, f, dfp);
     453        3598 :   t = gel(T,p);
     454        3598 :   if (typ(t) == t_INT)
     455             :   {
     456        2099 :     if (f == 1) return precomp(p, f, dfp);
     457         133 :     t = gel(T,p) = gclone(precomp(p, f, F));
     458             :   }
     459             :   /* t = precomp(p, 2, F) = HGMDATA2[p] */
     460        1632 :   if (f == 2)
     461         651 :     return (F == dfp)? t: Flv_red(t, upowuu(p, dfp));
     462             :   /* f = 1 */
     463         981 :   T = cgetg(p, t_VECSMALL); q = upowuu(p, dfp);
     464       10909 :   for (m = n = 1; m < p; m++, n += p+1) T[m] = uel(t, n) % q;
     465         981 :   return T;
     466             : }
     467             : 
     468             : /* W[N] = teich(N ^ (N * sign(VPOLGA[N]))) mod p^2 */
     469             : static GEN
     470       33936 : teichmodp2(GEN VPOLGA, ulong p, ulong p2)
     471             : {
     472             :   ulong pN;
     473             :   long c, l, N;
     474       33936 :   GEN W = cgetg_copy(VPOLGA,&l); /* W[N] = teich(N)^N */
     475      142875 :   for (N = 1, pN = p; N < l; N++, pN += p)
     476      108939 :     if ((c = VPOLGA[N]))
     477             :     {
     478       82255 :       long N0; (void)z_lvalrem(N, p, &N0);
     479       82255 :       if (c < 0) N0 = Fl_inv(N0 % p, p);
     480       82255 :       W[N] = Fl_powu(N0, pN, p2);
     481             :     }
     482       33936 :   return W;
     483             : }
     484             : /* GF[i] = i! mod p^2 */
     485             : static GEN
     486       33935 : get_GF(ulong p, ulong p2)
     487             : {
     488       33935 :   GEN F = cgetg(p, t_VECSMALL);
     489             :   ulong i;
     490    33631658 :   F[1] = 1; F[2] = 2; for (i = 3; i < p; i++) F[i] = Fl_mul(F[i-1], i, p2);
     491       36395 :   return F;
     492             : }
     493             : 
     494             : /* p odd, return v s.t. v[i] = 1/i mod p, i < p/2 */
     495             : static GEN
     496       33930 : Flv_inv_p2(ulong p)
     497             : {
     498       33930 :   long i, g = pgener_Fl(p), ph = (p >> 1);
     499       33934 :   GEN v = cgetg(ph + 1, t_VECSMALL);
     500       33935 :   pari_sp av = avma;
     501       33935 :   GEN w = cgetg(ph, t_VECSMALL);
     502       34607 :   w[1] = g;
     503    17050023 :   for (i = 2; i < ph; i++) w[i] = Fl_mul(w[i-1], g, p); /* w[i]=g^i */
     504    17384481 :   for (i = 1; i < ph; i++)
     505             :   {
     506    17350550 :     long x = w[i], y = w[ph - i]; /* g^(-i) = - g^(ph - i) */
     507    17350550 :     if (x > ph) v[p - x] = y; else v[x] = p - y; /* 1/(-x) = y, 1/x = -y */
     508             :   }
     509       33931 :   v[1] = 1; return gc_const(av, v);
     510             : }
     511             : 
     512             : /* H[i] = i * (H_i - ((p-1)! + 1) / p) mod p */
     513             : static GEN
     514       33926 : get_GH(GEN F, ulong p)
     515             : {
     516       33926 :   long i, g, ph = p >> 1;
     517       33926 :   GEN H = cgetg(p, t_VECSMALL), v = Flv_inv_p2(p);
     518       33936 :   H[p-1] = (F[p-1] + 1) / p; g = p - H[p-1];
     519    16244608 :   for (i = 1; i <= ph; i++)
     520             :   {
     521    16210683 :     long c = g = Fl_add(g, v[i], p);
     522    16214789 :     H[i] = c = Fl_mul(c, i, p);
     523    16213529 :     H[p-1 - i] = Fl_neg(Fl_add(c, g, p), p);
     524             :   }
     525       33925 :   return H;
     526             : }
     527             : /* p odd, GPB[m+1] = Gamma_p(m/(p-1)) mod p^2 */
     528             : static GEN
     529       33935 : doprecompmodp2(ulong p, ulong p2)
     530             : {
     531       33935 :   GEN F = get_GF(p, p2), H = get_GH(F, p), GPV = cgetg(p, t_VECSMALL);
     532             :   ulong r;
     533    32710966 :   for (r = 1; r < p; r++)
     534             :   {
     535    32677039 :     long t = Fl_mul(F[r], 1 + p * H[r], p2);
     536    32671386 :     GPV[p - r] = odd(r)? t: p2 - t;
     537             :   }
     538       33927 :   return GPV;
     539             : }
     540             : 
     541             : /***************************************************************/
     542             : /*        PART III: Motive initializations depending on p      */
     543             : /***************************************************************/
     544             : static GEN
     545       37597 : get_L1(GEN hgm, long PFM1, long f)
     546             : {
     547       37597 :   GEN VBEDEN = hgm_get_VBEDEN(hgm);
     548       37597 :   GEN VBENUM = hgm_get_VBENUM(hgm), v;
     549       37597 :   long DEG = hgm_get_DEG(hgm), j;
     550       37597 :   v = const_vecsmall(PFM1, f * hgm_get_TT(hgm));
     551      164024 :   for (j = 1; j <= DEG; j++)
     552      126427 :     if (PFM1 % VBEDEN[j] == 0)
     553             :     {
     554      118307 :       long t = ((-(PFM1 / VBEDEN[j]) * VBENUM[j]) % PFM1) + 1;
     555      118307 :       if (t <= 0) t += PFM1;
     556      118307 :       v[t] -= f;
     557             :     }
     558       37597 :   return v;
     559             : }
     560             : /* Stickelberger valuation: L0(p, f, m), power of (-p)
     561             :  * = f * OFFMPOL - sum_{e < f, N < l} gamma_N [N p^e m / (p^f-1)]*/
     562             : /* 1) amortized (all m < p^f-1), good when p >> N log N. Assume f=1 */
     563             : static GEN
     564       33963 : get_L0(GEN hgm, ulong PM1)
     565             : {
     566       33963 :   GEN perm, vt, vc, VL0, VPOLGA = hgm_get_VPOLGA(hgm);
     567       33962 :   long w, S, c, d, j, N, m, l = lg(VPOLGA), D = (l * (l-1)) >> 1;
     568       33962 :   vt = cgetg(D+1, t_VECSMALL);
     569       33962 :   vc = cgetg(D+1, t_VECSMALL);
     570      109087 :   for (N = 2, c = 1; N < l; N++)
     571       75126 :     if (VPOLGA[N])
     572             :     {
     573             :       ulong Q;
     574             :       long q;
     575      195489 :       for (q = 1, Q = 0; q <= N; q++, Q += PM1, c++)
     576             :       {
     577      147118 :         vt[c] = ceildivuu(Q, N);
     578      147117 :         vc[c] = VPOLGA[N];
     579             :       }
     580             :     }
     581       33961 :   setlg(vt,c); setlg(vc,c); perm = vecsmall_indexsort(vt);
     582       33968 :   vt = vecsmallpermute(vt, perm);
     583       33968 :   vc = vecsmallpermute(vc, perm);
     584       33969 :   w = vt[1];
     585      147157 :   for (j = 2, d = 1; j < c; j++)
     586      113188 :     if (vt[j] == w) vc[d] += vc[j];
     587             :     else
     588             :     {
     589       79301 :       d++;
     590       79301 :       vt[d] = w = vt[j];
     591       79301 :       vc[d] = vc[j];
     592             :     }
     593       33969 :   d++; vt[d] = PM1; vc[d] = 0; /* sentinels */
     594       33969 :   VL0 = cgetg(PM1+1, t_VECSMALL); S = hgm_get_OFFMPOL(hgm);
     595      147242 :   for (j = 1; j < d; j++)
     596             :   {
     597    34893408 :     for (m = vt[j]; m < vt[j+1]; m++) VL0[m+1] = S;
     598      113271 :     S -= vc[j+1];
     599             :   }
     600       33971 :   return VL0;
     601             : }
     602             : /* 2) direct computation (single m), good when p << N log N */
     603             : static long
     604      614870 : L0(GEN hgm, ulong p, long f, long PFM1, long m)
     605             : {
     606      614870 :   GEN VPOLGA = hgm_get_VPOLGA(hgm);
     607      614870 :   long e, S = hgm_get_OFFMPOL(hgm) * f, l = lg(VPOLGA);
     608             :   ulong pem;
     609     1858446 :   for (e = 0, pem = m; e < f; e++, pem *= p)
     610             :   {
     611     1243575 :     long N, s = 0;
     612             :     ulong Npem;
     613     9691550 :     for (N = 1, Npem = pem; N < l; N++, Npem += pem)
     614     8447975 :       if (VPOLGA[N]) s += VPOLGA[N] * (long)(Npem / PFM1);
     615     1243575 :     S -= s;
     616             :   }
     617      614871 :   return S;
     618             : }
     619             : 
     620             : static GEN
     621        2947 : get_teich1(GEN VPOLGA, GEN ZP, ulong p)
     622             : {
     623        2947 :   long l = lg(VPOLGA), N;
     624        2947 :   GEN v = zerovec(l - 1);
     625       12712 :   for (N = 2; N < l; N++)
     626        9765 :     if (VPOLGA[N])
     627             :     {
     628        3276 :       long N0; (void)z_lvalrem(N, p, &N0);
     629        3276 :       gel(v, N) = teich(gaddsg(N0, ZP));
     630             :     }
     631        2947 :   return v;
     632             : }
     633             : static GEN
     634        3661 : get_teich(GEN VPOLGA, GEN ZP, ulong p, long f, long PFM1)
     635             : {
     636             :   GEN v, Q;
     637             :   long l, N;
     638        3661 :   if (f == 1) return get_teich1(VPOLGA, ZP, p);
     639         714 :   l = lg(VPOLGA); v = zerovec(l - 1); Q = utoipos(PFM1 / (p - 1));
     640        4277 :   for (N = 2; N < l; N++)
     641        3563 :     if (VPOLGA[N])
     642             :     {
     643        1421 :       long N0; (void)z_lvalrem(N, p, &N0);
     644        1421 :       gel(v,N) = Qp_sqrtn(gdivsg(N0, teich(gaddsg(N0,ZP))), Q, NULL);
     645             :     }
     646         714 :   return v;
     647             : }
     648             : 
     649             : /* compute N^{-[N*s-1-(N*s-1)\p]}, N*s=a/(p^f-1) with a=N*m*p^e */
     650             : static GEN
     651      984447 : gapnpow(GEN T, long p, long f, long PFM1, long N, ulong Nmpe)
     652             : {
     653             :   GEN Vr;
     654             :   long Nm, i, S;
     655             : 
     656      984447 :   if (f == 1) return gpowgs(T, Nmpe);
     657      946295 :   Nm = Fl_neg(Nmpe, PFM1);
     658      946295 :   Vr = cgetg(f + 1, t_VECSMALL);
     659     3372390 :   for (i = 1; i <= f; i++) { Vr[i] = Nm % p; Nm /= p; }
     660     2426095 :   S = Vr[1]; for (i = 1; i < f; i++) S = Vr[f + 1 - i] + p * S;
     661      946295 :   return gdiv(gpowgs(T, S), powuu(N, Vr[1]));
     662             : }
     663             : 
     664             : /* prod_{j = 1}^DEG prod_{e = 0}^{f - 1}
     665             :  *   gamma_p({p^e * (VAL[j] + m/(p^f-1))}) /
     666             :  *   gamma_p({p^e * (VBE[j] + m/(p^f-1))}), a p-unit. 0 <= m < p^f-1 */
     667             : static GEN
     668      235082 : hgmC(GEN VPOLGA, GEN GPV, GEN TEICH, ulong p, long f, ulong PFM1, ulong m, long dfp)
     669             : {
     670      235082 :   GEN r = gen_1;
     671      235082 :   long c, e, N, l = lg(VPOLGA);
     672             :   ulong Nmpe, mpe;
     673      700169 :   for (e = 0, mpe = m; e < f; e++, mpe = Fl_mul(mpe, p, PFM1))
     674     3819412 :     for (N = 1, Nmpe = mpe; N < l; N++, Nmpe = Fl_add(Nmpe, mpe, PFM1))
     675     3354322 :       if ((c = VPOLGA[N]))
     676             :       { /* Gamma_p(frac(N*m*p^e/(p^f-1)))  */
     677     1449547 :         GEN z = utoi(GPV[Nmpe + 1]);
     678     1449551 :         if (N > 1) z = gmul(z, gapnpow(gel(TEICH, N), p, f, PFM1, N, Nmpe));
     679             :         /* z = prod_{0 <= j < N} Gamma_p( frac(p^e * (m/(p^f-1) + j/N)) ) */
     680     1449536 :         r = gmul(r, gpowgs(z, c));
     681             :       }
     682      235062 :   if (typ(r) != t_PADIC) r = cvtop(r, utoipos(p), dfp);
     683      235060 :   return r;
     684             : }
     685             : 
     686             : /* Same modulo p^2, Wm[N] = teich(N^(sign(VPOLGA[N]) * m)) */
     687             : static ulong
     688    16369319 : hgmCmodp2(GEN VPOLGA, GEN Wm, GEN GPV, ulong PFM1, ulong m, ulong p2)
     689             : {
     690    16369319 :   long l = lg(VPOLGA), c, N;
     691    16369319 :   ulong res = 1, Nm;
     692    64920597 :   for (N = 1, Nm = m; N < l; N++, Nm = Fl_add(Nm, m, PFM1))
     693    48497818 :     if ((c = VPOLGA[N]))
     694             :     {
     695             :       ulong z;
     696    37049729 :       if (c > 0) z = GPV[Nm + 1];
     697             :       else
     698             :       {
     699    18994951 :         c = -c;
     700    18994951 :         if (!Nm) z = 1;
     701             :         else
     702             :         {
     703    18975548 :           z = GPV[PFM1 + 1 - Nm];
     704    18975548 :           if (!odd(Nm)) z = p2 - z; /* GPV[Nm+1]^(-1) by reflection formula */
     705             :         }
     706             :       }
     707    37716468 :       if (N > 1) z = Fl_mul(z, Wm[N], p2);
     708    37690844 :       z = Fl_powu(z, c, p2);
     709             :       /* z = (GPV[Nm+1] * teich(N^m))^VPOLGA[N] */
     710    37291570 :       res = res == 1? z: Fl_mul(res, z, p2);
     711             :     }
     712    16358035 :   return res;
     713             : }
     714             : 
     715             : /***************************************************************/
     716             : /*        PART IV: Motive functions depending on p, t          */
     717             : /***************************************************************/
     718             : static long
     719       37891 : get_dfp(GEN hgm, long p, long f)
     720             : {
     721       37891 :   long DEG = hgm_get_DEG(hgm), WT = hgm_get_WT(hgm);
     722       37889 :   long dfp = ceil(log((double)2*DEG) / log((double)p) + f * (WT * 0.5 ) + 1e-5);
     723       37889 :   if (p == 2) dfp++; /* FIXME: why ? */
     724       37889 :   return dfp;
     725             : }
     726             : /* V being a vecsmall, compute all p^TT*hgmC(m)/hgmC(0) for indices in V */
     727             : static GEN
     728        3668 : hgmCall(GEN hgm, ulong p, long f, long dfp, GEN V)
     729             : {
     730        3668 :   GEN v, c, GPV, VL0, VL1, TEICH, ZP = zeropadic_shallow(utoipos(p), dfp);
     731        3668 :   GEN VPOLGA = hgm_get_VPOLGA(hgm);
     732        3668 :   ulong i, PFM1 = upowuu(p, f) - 1, lV = V? lg(V): PFM1+1;
     733        3668 :   long l0, fTT = f * hgm_get_TT(hgm);
     734             : 
     735        3668 :   GPV = doprecomp(p, f, dfp);
     736        3661 :   VL0 = (V && f == 1)? get_L0(hgm, PFM1): NULL;
     737        3661 :   VL1 = get_L1(hgm, PFM1, f);
     738        3661 :   TEICH = get_teich(VPOLGA, ZP, p, f, PFM1);
     739        3661 :   l0 = hgm_get_OFFMPOL(hgm) * f;
     740        3661 :   if (V)
     741             :   {
     742          42 :     v = cgetg(lV, t_VEC);
     743          42 :     i = 1;
     744             :   }
     745             :   else
     746             :   {
     747        3619 :     v = cgetg(lV+1, t_POL); v[1] = evalsigne(1)|evalvarn(0);
     748        3619 :     gel(v,2) = c = powuu(p, fTT); /* m = 0 */
     749        3619 :     i = 2;
     750             :   }
     751      618606 :   for (; i < lV; i++)
     752             :   {
     753      614947 :     long m = V? V[i]: i-1;
     754      614947 :     long L = VL0? VL0[m+1]: L0(hgm, p, f, PFM1, m);
     755      614948 :     long e = L + VL1[m+1];
     756      614948 :     if (!V && e >= dfp) c = gen_0;
     757             :     else
     758             :     { /* FIXME: dfp is fishy in Jordantame, don't trust if V = NULL */
     759      235079 :       pari_sp av = avma;
     760      235079 :       c = hgmC(VPOLGA, GPV, TEICH, p, f, PFM1, m, dfp);
     761      235059 :       setvalp(c, e); if (odd(L ^ l0)) c = gneg(c);
     762      235058 :       c = gerepileupto(av, padic_to_Q(c));
     763             :     }
     764      614945 :     gel(v, V? i: i+1) = c;
     765             :   }
     766        3659 :   return v;
     767             : }
     768             : /* Same mod p^2, f = 1 */
     769             : static GEN
     770       33935 : hgmCallmodp2(GEN hgm, ulong p)
     771             : {
     772       33935 :   GEN C, GPV, VL0, VL1, W1, Wm, VPOLGA = hgm_get_VPOLGA(hgm);
     773       33935 :   long l = lg(VPOLGA), TT = hgm_get_TT(hgm);
     774       33935 :   ulong m, PFM1 = p - 1, p2 = p * p;
     775             : 
     776       33935 :   if (p & HIGHMASK) pari_err_OVERFLOW("hgmCallmodp2");
     777       33928 :   VL0 = get_L0(hgm, PFM1);
     778       33936 :   VL1 = get_L1(hgm, PFM1, 1);
     779       33936 :   W1 = teichmodp2(VPOLGA, p, p2);
     780       33936 :   Wm = shallowcopy(W1);
     781       33935 :   GPV = doprecompmodp2(p, p2);
     782       33927 :   C = cgetg(PFM1+2, t_VECSMALL); C[1] = evalvarn(0);
     783       33928 :   C[2] = TT > 1? 0: (TT == 1? p : 1); /* m = 0 */
     784    32080578 :   for (m = 1; m < PFM1; m++)
     785             :   {
     786    32547150 :     long e = VL0[m + 1] + VL1[m + 1], j;
     787             :     ulong c;
     788    32547150 :     if (e >= 2) c = 0;
     789             :     else
     790             :     {
     791    16315354 :       c = hgmCmodp2(VPOLGA, Wm, GPV, PFM1, m, p2);
     792    16349773 :       if (odd(VL0[m + 1] ^ VL0[1])) c = Fl_neg(c, p2);
     793    15843370 :       if (e == 1) c = (c % p) * p;
     794             :     }
     795    32075166 :     C[m + 2] = c;
     796    95220479 :     for (j = 2; j < l; j++)
     797    63173829 :       if (VPOLGA[j]) Wm[j] = Fl_mul(Wm[j], W1[j], p2);
     798             :   }
     799       29088 :   return C;
     800             : }
     801             : 
     802             : /* 1 / (1-q) ~ 1 + q + ... + q^n; n >= 1 */
     803             : static ulong
     804        3906 : inv(ulong q, long n)
     805             : {
     806        3906 :   ulong z = q + 1;
     807             :   long i;
     808        8064 :   for (i = 1; i < n; i++) z = z * q + 1;
     809        3906 :   return z;
     810             : }
     811             : 
     812             : /* General H function: C(teich(f + O(p^dfp))) / (1 - p^f) */
     813             : static GEN
     814        3906 : hgmH(GEN C, long p, long f, long dfp, GEN t)
     815             : {
     816        3906 :   GEN q = powuu(p, dfp), z;
     817             :   long n;
     818        3906 :   (void)Q_lvalrem(t, p, &t);
     819        3906 :   z = Rg_to_Fp(t, q);
     820        3906 :   z = Zp_teichmuller(z, utoipos(p), dfp, q);
     821        3906 :   z = FpX_eval(C, z, q);
     822        3906 :   n = dfp / f; if (!(dfp % f)) n--;
     823        3906 :   z = Fp_mulu(z, inv(upowuu(p,f), n), q);
     824        3906 :   return Fp_center(z, q,  shifti(q,-1));
     825             : }
     826             : static GEN
     827       33927 : hgmHmodp2(GEN C, ulong p, GEN t)
     828             : {
     829       33927 :   ulong p2 = p * p, wt = Fl_powu(Rg_to_Fl(t, p2), p, p2);
     830       33934 :   return stoi( Fl_center(Fl_mul(Flx_eval(C, wt, p2), 1+p, p2), p2, p2 >> 1) );
     831             : }
     832             : 
     833             : enum { C_OK = 0, C_FAKE, C_BAD, C_TAME0, C_TAME1};
     834             : static GEN hgmU(GEN hgm, long p, long f, GEN t, long dfp);
     835             : static long hgmclass(GEN hgm, long p, GEN t);
     836             : static GEN
     837       37849 : hgmtrace(GEN hgm, long p, long f, GEN t, long c)
     838             : {
     839       37849 :   long dfp = get_dfp(hgm, p, f);
     840       37847 :   if (c == C_FAKE) return hgmU(hgm, p, f, t, dfp);
     841       37560 :   if (f == 1 && dfp <= 2) return hgmHmodp2(hgmCallmodp2(hgm, p), p, t);
     842        3625 :   return hgmH(hgmCall(hgm, p, f, dfp, NULL), p, f, dfp, t);
     843             : }
     844             : 
     845             : static GEN
     846         700 : F2v_factorback(GEN E)
     847             : {
     848         700 :   long i, l = lg(E);
     849         700 :   GEN c = gen_1;
     850        4235 :   for (i = 1; i < l; i++)
     851        3535 :     if (odd(E[i])) c = muliu(c, i);
     852         700 :   return c;
     853             : }
     854             : 
     855             : static long
     856       28419 : Q_krois(GEN T, long p) { return krouu(Rg_to_Fl(T, (p == 2)? 8: p), p); }
     857             : static long
     858       33931 : hgmsign(GEN hgm, long p, GEN t)
     859             : {
     860             :   GEN u;
     861       33931 :   if (odd(hgm_get_WT(hgm))) return 1;
     862       28418 :   t = ginv(t); u = gmul(gsubsg(1, t), hgm_get_SIGNPAR(hgm));
     863       28419 :   return odd(hgm_get_DEG(hgm))? -Q_krois(u, p): Q_krois(gmul(gneg(t), u), p);
     864             : }
     865             : /* conductor of the central character */
     866             : static GEN
     867         259 : hgmcharcond(GEN hgm, GEN t)
     868             : {
     869             :   GEN u;
     870         259 :   if (odd(hgm_get_WT(hgm))) return gen_1;
     871         210 :   t = ginv(t); u = gmul(gsubsg(1, t), hgm_get_SIGNPAR(hgm));
     872         210 :   if (!odd(hgm_get_DEG(hgm))) u = gmul(gneg(t), u);
     873         210 :   if (typ(u) == t_FRAC) u = mulii(gel(u,1), gel(u,2));
     874         210 :   return gequal0(u) ? gen_1 : coredisc(u);
     875             : }
     876             : /* valuations of central character conductor at BAD */
     877             : static GEN
     878          98 : get_achi(GEN hgm, GEN t, GEN BAD)
     879             : {
     880          98 :   long i, l = lg(BAD);
     881          98 :   GEN Nchi = hgmcharcond(hgm, t), v = cgetg(l, t_VECSMALL);
     882         231 :   for (i = 1; i < l; i++) v[i] = Z_lval(Nchi, BAD[i]);
     883          98 :   return v;
     884             : }
     885             : 
     886             : /* [a,a*r, ..., a*r^n] */
     887             : static GEN
     888         420 : upowers_u(ulong r, long n, ulong a)
     889             : {
     890         420 :   long i, l = n + 2, t = a;
     891         420 :   GEN v = cgetg(l, t_VECSMALL);
     892        2562 :   for (i = 1; i < l; i++) { v[i] = t; t *= r; }
     893         420 :   return v;
     894             : }
     895             : static GEN
     896       36873 : powers_u(ulong r, long n, ulong a)
     897             : {
     898       36873 :   long i, l = n + 2;
     899       36873 :   GEN v = cgetg(l, t_VEC), t = utoi(a);
     900      106585 :   for (i = 1; i < l; i++) { gel(v,i) = t; t = muliu(t, r); }
     901       36874 :   return v;
     902             : }
     903             : static GEN
     904       36872 : mkpowers(long p, long d, long WT)
     905             : {
     906             :   ulong q, r;
     907       36872 :   if (WT == 0) q = r = 1;
     908       35696 :   else if (!odd(d)) q = r = upowuu(p, WT);
     909       21888 :   else if (WT == 1) { q = 1; r = p; }
     910       21888 :   else { q = upowuu(p, WT >> 1); r = q*q; if (odd(WT)) r *= p; }
     911       36872 :   return powers_u(r, (d-1)>>1, q);
     912             : }
     913             : 
     914             : /* complete local factor E (of degree d) at p using functional equation
     915             :  * E(T) = SIGN p^(WT*d)/2 T^d E(1/(p^WT*T)) */
     916             : static GEN
     917       36871 : Efuneq(GEN E, long p, long d, long WT, long SIGN, long B)
     918             : {
     919       36871 :   long j = maxss(d - B, 0), k = d + 1 - j, nE = lg(E)-1, l = (d + 1) >> 1;
     920       36871 :   GEN T = cgetg(k + 1, t_VEC), vp = mkpowers(p, d, WT);
     921       39708 :   for (; j < l; j++, k--)
     922             :   {
     923        2835 :     GEN q = gel(vp,l-j);
     924        2835 :     if (SIGN < 0) q = negi(q); /* q = SIGN * p^(WT(d-2*j) / 2) */
     925        2835 :     gel(T, k) = gmul(q, RgX_coeff(E, j));
     926             :   }
     927       39497 :   for (; k >= nE; k--) gel(T, k) = gen_0;
     928      108604 :   for (; k > 0; k--) gel(T, k) = gel(E, k+1);
     929       36873 :   return RgV_to_RgX(T,0);
     930             : }
     931             : 
     932             : static long
     933       37641 : hgmclass(GEN hgm, long p, GEN t)
     934             : {
     935       37641 :   GEN BAD = hgm_get_BAD(hgm);
     936             :   long ap, bp;
     937       37640 :   if (!umodiu(BAD, p))
     938             :   {
     939         679 :     long v = Q_lval(t, p);
     940         679 :     if (v && v + Q_lval(hgm_get_MVALUE(hgm), p) == 0)
     941             :     {
     942         161 :       GEN N = hgmcharcond(hgm, t);
     943         161 :       if (umodiu(N, p)) return C_FAKE; /* a(chi) = 0 */
     944             :     }
     945         553 :     return C_BAD;
     946             :   }
     947       36961 :   if (typ(t) == t_INT)
     948             :   {
     949       23993 :     ap = umodiu(t, p); if (!ap) return C_TAME0;
     950       23965 :     bp = 1;
     951             :   }
     952             :   else
     953             :   {
     954       12968 :     ap = umodiu(gel(t,1), p); if (!ap) return C_TAME0;
     955       12961 :     bp = umodiu(gel(t,2), p); if (!bp) return C_TAME0;
     956             :   }
     957       36863 :   return ap == bp? C_TAME1 : C_OK;
     958             : }
     959             : 
     960             : /* p good or Tame1; return local factor at p: 1/E + O(x^(B+1)); t a t_VEC,
     961             :  * C a t_VECSMALL giving their class. */
     962             : static GEN
     963       36962 : frobpoltrunc(GEN hgm, GEN t, long c, long p, long B, long *pF)
     964             : {
     965       36962 :   long DEG = hgm_get_DEG(hgm), WT = hgm_get_WT(hgm);
     966       36961 :   long D = isint1(t)? (odd(WT) ? DEG - 2 : DEG - 1): DEG;
     967       36960 :   long f, mi, m, q = upowuu(p, WT >> 1);
     968             :   GEN E, s;
     969             : 
     970       36960 :   mi = minss(B, (c == C_FAKE)? D: D >> 1);
     971       36960 :   m = (mi == D && c == C_TAME1 && !odd(DEG))? mi: mi+1;
     972       36960 :   s = cgetg(m+1, t_POL); s[1] = evalsigne(1)|evalvarn(0);
     973       74797 :   for (f = 1; f < m; f++) gel(s, f+1) = negi( hgmtrace(hgm, p, f, t, c) );
     974       36952 :   s = RgX_renormalize_lg(s, m+1);
     975       36951 :   *pF = 0; s = RgXn_expint(s, m);
     976       36948 :   if (mi == D) return s;
     977       36871 :   if (c == C_TAME1)
     978             :   {
     979        2940 :     long SIGN = kroiu(hgm_get_U(hgm), p);
     980        2939 :     if (odd(DEG))
     981         203 :       E = Efuneq(s, p, DEG-1, WT, SIGN, B);
     982             :     else
     983             :     {
     984        2736 :       GEN T = deg1pol_shallow(stoi(- SIGN * q), gen_1, 0);
     985        2737 :       E = RgXn_mul(s, RgXn_inv(T, m), m);
     986        2737 :       E = Efuneq(E, p, DEG - 2, WT, 1, B);
     987        2737 :       if (!gequal1(t) || !odd(WT)) E = gmul(E, T);
     988             :     }
     989        2940 :     *pF = 1;
     990        2940 :     if (!odd(WT))
     991             :     {
     992             :       GEN T, u, t0;
     993         217 :       long v = Q_lvalrem(gsubgs(t, 1), p, &t0);
     994         217 :       if (!odd(v))
     995             :       {
     996          63 :         if (typ(t0) == t_FRAC) t0 = mulii(gel(t0,1), gel(t0,2));
     997          63 :         u = coredisc(mulii(t0, hgm_get_U0(hgm)));
     998          63 :         T = deg1pol_shallow(stoi(-kroiu(u, p) * q), gen_1, 0);
     999          63 :         E = gmul(E, T); *pF = 0;
    1000             :       }
    1001             :     }
    1002             :   }
    1003             :   else
    1004             :   {
    1005       33931 :     long SIGN = hgmsign(hgm, p, t);
    1006       33931 :     E = Efuneq(s, p, DEG, WT, SIGN, B);
    1007             :   }
    1008       36872 :   return E;
    1009             : }
    1010             : 
    1011             : GEN
    1012          84 : hgmcoef(GEN hgm, GEN t, GEN n)
    1013             : {
    1014          84 :   pari_sp av = avma;
    1015          84 :   GEN T, P, E, F = check_arith_all(n, "hgmcoef");
    1016             :   long i, lP;
    1017             : 
    1018          84 :   if (!checkhgm(hgm)) pari_err_TYPE("hgmcoef", hgm);
    1019          84 :   if (!is_rational_t(typ(t))) pari_err_TYPE("hgmcoef",t);
    1020          77 :   if (hgm_get_SWAP(hgm)) t = ginv(t);
    1021          77 :   if (!F) { F = Z_factor(n); P = gel(F,1); }
    1022             :   else
    1023             :   {
    1024          21 :     P = gel(F,1);
    1025          21 :     if (lg(P) == 1 || signe(gel(P,1)) <= 0) return gen_0;
    1026          21 :     n = typ(n) == t_VEC? gel(n,1): factorback(F);
    1027             :   }
    1028          77 :   if (signe(n) <= 0) pari_err_DOMAIN("hgmcoef", "n", "<=", gen_0, n);
    1029          77 :   E = gel(F,2); lP = lg(P); T = gen_1;
    1030          77 :   T = gen_1;
    1031         161 :   for (i = 1; i < lP; i++)
    1032             :   {
    1033          84 :     long e, p = itos(gel(P, i)), f = itos(gel(E, i)), c = hgmclass(hgm, p, t);
    1034             :     GEN A;
    1035          84 :     if (c == C_BAD) pari_err_IMPL("hgmcoef for bad primes");
    1036          84 :     A = frobpoltrunc(hgm, t, c, p, f, &e);
    1037          84 :     T = gmul(T, RgX_coeff(RgXn_inv(A, f+1), f));
    1038             :   }
    1039          77 :   return gerepilecopy(av, T);
    1040             : }
    1041             : 
    1042             : static GEN
    1043          14 : count2list(GEN E)
    1044             : {
    1045          14 :   long i, j, r, l = lg(E);
    1046          14 :   GEN v = cgetg(zv_sum(E)+1, t_VECSMALL);
    1047          49 :   for (i = j = 1; i < l; i++)
    1048          56 :     if ((r = E[i])) while(r--) v[j++] = i;
    1049          14 :   setlg(v,j); return v;
    1050             : }
    1051             : static GEN hgminit_i(GEN val, GEN vbe);
    1052             : #if 0
    1053             : static long
    1054             : minval(GEN F, long p)
    1055             : {
    1056             :   long i, d = degpol(F), m = LONG_MAX;
    1057             :   for (i = 1; i <= d; i++) m = minss(m, Z_lval(gel(F,i+2), p) / i);
    1058             :   return m;
    1059             : }
    1060             : static GEN
    1061             : cycloE_filter(GEN A, long p)
    1062             : {
    1063             :   long i, j, l = lg(A);
    1064             :   GEN B = shallowcopy(A);
    1065             :   if (p < l)
    1066             :     for (i = p; i < l; i += p)
    1067             :       if (A[i])
    1068             :       {
    1069             :         (void)z_lvalrem(i / p, p, &j);
    1070             :         B[j] += A[i];
    1071             :         B[i] = 0;
    1072             :       }
    1073             :   return B;
    1074             : }
    1075             : /* doesn't work because A/B is not valid hypergeometric data */
    1076             : static GEN
    1077             : eulfacbadnew(GEN hgm, GEN t, long p, long *pe)
    1078             : {
    1079             :   GEN AB = hgm_get_CYCLOE(hgm), gp = utoipos(p);
    1080             :   GEN A = cycloE_filter(gel(AB,1), p);
    1081             :   GEN B = cycloE_filter(gel(AB,2), p);
    1082             :   GEN hgmp = hgminit_i(count2list(A), count2list(B));
    1083             :   GEN cF, E, F = hgmeulerfactor(hgmp, t, p, &E);
    1084             :   long s, v, d = degpol(F), w = hgm_get_WT(hgm);
    1085             :   *pe = 0;
    1086             :   if (!d) return pol_1(0);
    1087             :   cF = leading_coeff(F);
    1088             :   v = Z_lvalrem(cF, p, &cF);
    1089             :   if (!is_pm1(cF)) return pol_1(0);
    1090             :   s = minval(F, p);
    1091             :   if (s) { F = RgX_unscale(F, powis(gp, -s)); v -= s * d; }
    1092             :   if ((2 * v) % d || (!odd(w) && v % d)) return F;
    1093             :   return RgX_unscale(F, powis(gp, w / 2 - v / d));
    1094             : }
    1095             : #endif
    1096             : static GEN Jordantameexpo(GEN hgm, long v, GEN t0, long p, long *pe);
    1097             : static GEN
    1098       37271 : hgmeulerfactorlimit(GEN hgm, GEN t, long p, long d, long flag, long *pe)
    1099             : {
    1100       37271 :   long c = hgmclass(hgm, p, t);
    1101       37269 :   if (c == C_TAME0)
    1102             :   {
    1103          91 :     long v = Q_lvalrem(t, p, &t);
    1104          91 :     return Jordantameexpo(hgm, v, t, p, pe);
    1105             :   }
    1106       37178 :   if (c != C_BAD) return frobpoltrunc(hgm, t, c, p, d, pe);
    1107         301 :   if (flag) { *pe = -1; return gen_0; } else { *pe = 0; return pol_1(0); }
    1108             : }
    1109             : 
    1110             : GEN
    1111         455 : hgmeulerfactor(GEN hgm, GEN t, long p, GEN* pE)
    1112             : {
    1113         455 :   pari_sp av = avma;
    1114             :   long e, B;
    1115             :   GEN P;
    1116         455 :   if (!checkhgm(hgm)) pari_err_TYPE("hgmeulerfactor", hgm);
    1117         455 :   if (!is_rational_t(typ(t))) pari_err_TYPE("hgmeulerfactor",t);
    1118         448 :   if (hgm_get_SWAP(hgm)) t = ginv(t);
    1119         448 :   B = (long)(hgm_get_DEG(hgm) * log(p)) + 1;
    1120         448 :   P = gerepilecopy(av, hgmeulerfactorlimit(hgm, t, p, B, 1, &e));
    1121         434 :   if (pE) *pE = stoi(e);
    1122         434 :   return P;
    1123             : }
    1124             : 
    1125             : /***********************************************************************/
    1126             : /*                        Tame Euler factors                           */
    1127             : /***********************************************************************/
    1128             : /* FIXME: implement properly like RgXn_sqrt */
    1129             : static GEN
    1130          42 : RgXn_sqrtnu(GEN h, long f, long e)
    1131             : {
    1132          42 :   if (f == 1) return h;
    1133           7 :   if (f == 2) return RgXn_sqrt(h, e);
    1134           0 :   return ser2rfrac_i(gsqrtn(RgX_to_ser(h, e + 2), utoipos(f), NULL, 0));
    1135             : }
    1136             : static GEN
    1137          77 : Jordantame(GEN hgm, GEN t0, long m, long p)
    1138             : {
    1139             :   GEN P, T, C, ZP, V;
    1140             :   long d, phim, f, j, c, q, qm, dfp;
    1141             : 
    1142          77 :   if (m == 1)
    1143             :   {
    1144          35 :     long e = hgm_get_WT(hgm) - get_b1(hgm_get_CYCLOE(hgm));
    1145          35 :     return deg1pol_shallow(negi(powuu(p, (e+1) >> 1)), gen_1, 0);
    1146             :   }
    1147          42 :   phim = eulerphiu(m); f = Fl_order(p % m, phim, m);
    1148          42 :   d = phim + 1; /* DEG >= phim >= f */
    1149          42 :   q = upowuu(p, f); qm = (q - 1) / m;
    1150          42 :   V = cgetg(phim + 1, t_VECSMALL);
    1151         133 :   for (j = c = 1; j < m; j++)
    1152          91 :     if (cgcd(j, m) == 1) V[c++] = j * qm;
    1153          42 :   dfp = get_dfp(hgm, p, f);
    1154          42 :   C = hgmCall(hgm, p, f, dfp, V);
    1155          42 :   ZP = zeropadic_shallow(utoipos(p), dfp);
    1156          42 :   T = teich(gadd(t0, ZP)); P = pol_1(0);
    1157         133 :   for (j = 1; j < lg(V); j++)
    1158             :   {
    1159          91 :     GEN Q, c = gmul(gpowgs(T, V[j]), gel(C,j));
    1160          91 :     Q = RgXn_red_shallow(RgX_shift_shallow(RgX_Rg_mul(P, c), f), d);
    1161          91 :     P = RgX_sub(P, Q); /* P * (1 - c x^f) mod X^d */
    1162             :   }
    1163          42 :   return centerlift(RgXn_sqrtnu(P, f, d));
    1164             : }
    1165             : 
    1166             : static GEN
    1167         105 : eulfactameinit(GEN hgm, long v)
    1168             : {
    1169         105 :   GEN C = hgm_get_CYCLOE(hgm);
    1170         105 :   if (!v) pari_err_BUG("hgmeulerfactor [incorrect t in eulfactame]");
    1171         105 :   if (hgm_get_SWAP(hgm)) v = -v;
    1172         105 :   return gel(C, (v < 0)? 2: 1);
    1173             : }
    1174             : static long
    1175          14 : tameexpo(GEN hgm, long v)
    1176             : {
    1177          14 :   GEN W = eulfactameinit(hgm,v);
    1178          14 :   long e, m, l = lg(W);
    1179          42 :   for (m = 1, e = 0; m < l; m++)
    1180          28 :     if (W[m] && v % m == 0) e += eulerphiu(m);
    1181          14 :   return hgm_get_DEG(hgm) - e;
    1182             : }
    1183             : static GEN
    1184          91 : Jordantameexpo(GEN hgm, long v, GEN t0, long p, long *pe)
    1185             : {
    1186          91 :   GEN P = pol_1(0), W = eulfactameinit(hgm, v);
    1187          91 :   long e, m, l = lg(W);
    1188         308 :   for (m = 1, e = 0; m < l; m++)
    1189         217 :     if (W[m] && v % m == 0)
    1190             :     {
    1191          77 :       P = gmul(P, Jordantame(hgm, t0, m, p));
    1192          77 :       e += eulerphiu(m);
    1193             :     }
    1194          91 :   *pe = hgm_get_DEG(hgm) - e; return P;
    1195             : }
    1196             : 
    1197             : /***************************************************************/
    1198             : /*         PART IV.5: Fake wild primes for t                   */
    1199             : /***************************************************************/
    1200             : /* Compute g_q(r)=pi^sq(r)*gauss(r) */
    1201             : static GEN
    1202       57932 : gkgauss(long p, long f, GEN vp, long r, GEN ZP, GEN *sp)
    1203             : {
    1204       57932 :   GEN S = gen_0, P = gen_m1;
    1205       57932 :   long i, qm1 = vp[f+1] - 1;
    1206      272062 :   for (i = 1; i <= f; i++)
    1207             :   {
    1208      214130 :     GEN t = gfrac(sstoQ(r * vp[i], qm1));
    1209      214130 :     S = gadd(S, t);
    1210      214130 :     P = gmul(P, Qp_gamma(gadd(t, ZP)));
    1211             :   }
    1212       57932 :   *sp = gmulsg(p - 1, S); /* t_INT */
    1213       57932 :   return P;
    1214             : }
    1215             : 
    1216             : static GEN
    1217       12292 : hgmG(GEN hgm, long p, long f, GEN vp, long r, GEN ZP)
    1218             : {
    1219       12292 :   GEN S = gen_0, P = gen_1, VPOLGA = hgm_get_VPOLGA(hgm);
    1220             :   long n, c;
    1221       81088 :   for (n = 1; n < lg(VPOLGA); n++)
    1222       68796 :     if ((c = VPOLGA[n]))
    1223             :     {
    1224       57932 :       GEN sq, g = gkgauss(p, f, vp, r*n, ZP, &sq);
    1225       57932 :       S = addii(S, mulsi(c, sq));
    1226       57932 :       P = gmul(P, gpowgs(g, c));
    1227             :     }
    1228             :   /* p - 1 | S */
    1229       12292 :   return gmul(P, powis(utoineg(p), itos(S) / (p - 1)));
    1230             : }
    1231             : 
    1232             : /* multiplicity of -r / (q-1) in beta */
    1233             : static long
    1234       24584 : hgmmulti(GEN B, long q, long r)
    1235             : {
    1236       24584 :   long d = (q-1) / ugcd(r, q-1);
    1237       24584 :   return d >= lg(B)? 0: B[d];
    1238             : }
    1239             : 
    1240             : /* We must have w(M)^r*hgmQ(q,r)=hgmC(q,r)/hgmC(q,0) */
    1241             : static GEN
    1242       12292 : hgmQ(GEN hgm, long p, long f, GEN vp, long r, GEN ZP)
    1243             : {
    1244       12292 :   pari_sp av = avma;
    1245       12292 :   GEN B = gel(hgm_get_CYCLOE(hgm), 2);
    1246       12292 :   long q = vp[f+1], m0 = hgmmulti(B, q, 0), m1 = hgmmulti(B, q, r);
    1247       12292 :   GEN c = powis(utoipos(q), hgm_get_TT(hgm) + m0 - m1);
    1248       12292 :   if (odd(m0)) c = negi(c);
    1249       12292 :   return gerepileupto(av, padic_to_Q(gmul(c, hgmG(hgm, p, f, vp, r, ZP))));
    1250             : }
    1251             : 
    1252             : static GEN
    1253         287 : hgmU(GEN hgm, long p, long f, GEN t, long dfp)
    1254             : {
    1255         287 :   pari_sp av = avma;
    1256         287 :   GEN ZP = zeropadic_shallow(utoipos(p), dfp), vp = upowers_u(p, f, 1);
    1257         287 :   long q = vp[f+1], i;
    1258         287 :   GEN Q = cgetg(q+1, t_POL);
    1259         287 :   Q[1] = evalsigne(1)|evalvarn(0);
    1260       12579 :   for (i = 2; i <= q; i++) gel(Q, i) = hgmQ(hgm, p, f, vp, i-2, ZP);
    1261         287 :   t = p == 2? gen_1: gmul(hgm_get_MVALUE(hgm), t);
    1262         287 :   return gerepileupto(av, hgmH(Q, p, f, dfp, t));
    1263             : }
    1264             : 
    1265             : /***************************************************************/
    1266             : /*         PART V: Utility programs and main init              */
    1267             : /***************************************************************/
    1268             : static GEN
    1269           7 : cycloE2cyclo(GEN A, GEN B)
    1270           7 : { return mkvec2(count2list(A), count2list(B)); }
    1271             : #if 0
    1272             : GEN
    1273             : hgmalphatocyclo(GEN val, GEN vbe)
    1274             : { GEN C = get_CYCLOE(val,vbe); return cycloE2cyclo(gel(C,1), gel(C,2)); }
    1275             : #endif
    1276             : 
    1277             : static GEN
    1278      843850 : allprims(long n, GEN cache)
    1279             : {
    1280             :   long l, i, c;
    1281             :   GEN z, v;
    1282      843850 :   if (gel(cache,n)) return gel(cache,n);
    1283      418635 :   z = coprimes_zv(n); l = lg(z); v = cgetg(l, t_VEC);
    1284     3036782 :   for (i = c = 1; i < l; i++)
    1285     2618147 :     if (z[i]) gel(v, c++) = mkfracss(i, n);
    1286      418635 :   setlg(v, c); return gel(cache,n) = v;
    1287             : }
    1288             : static GEN
    1289      207956 : zv_to_prims(GEN D, GEN cache)
    1290             : {
    1291      207956 :   long i, l = lg(D);
    1292      207956 :   GEN v = cgetg(l, t_VEC);
    1293     1051806 :   for (i = 1; i < l; i++)
    1294             :   {
    1295      843857 :     if (D[i] <= 0) pari_err_TYPE("hgmcyclotoalpha", D);
    1296      843850 :     gel(v, i) = allprims(D[i], cache);
    1297             :   }
    1298      207949 :   return shallowconcat1(v);
    1299             : }
    1300             : static void
    1301      103978 : hgmcyclotoalpha(GEN *pA, GEN *pB)
    1302             : {
    1303      103978 :   GEN v, D = *pA, E = *pB;
    1304      103978 :   if (typ(D) != t_VECSMALL) D = gtovecsmall(D);
    1305      103978 :   if (typ(E) != t_VECSMALL) E = gtovecsmall(E);
    1306      103978 :   v = const_vec(maxss(vecsmall_max(D), vecsmall_max(E)), NULL);
    1307      103978 :   gel(v,1) = mkvec(gen_0);
    1308      103978 :   *pA = zv_to_prims(D, v);
    1309      103978 :   *pB = zv_to_prims(E, v);
    1310      103971 :   if (lg(*pA) != lg(*pB))
    1311           0 :     pari_err_TYPE("hgminit [incorrect lengths]", mkvec2(D,E));
    1312      103971 : }
    1313             : 
    1314             : static GEN
    1315         231 : hgmalphatogamma(GEN val, GEN vbe) { return get_VPOLGA(get_CYCLOE(val, vbe)); }
    1316             : 
    1317             : static void
    1318          21 : hgmgammatocyclo(GEN VPOLGA, GEN *pD, GEN *pE)
    1319             : {
    1320          21 :   long i, cn, cd, l = lg(VPOLGA);
    1321          21 :   GEN d, n, v = zero_zv(l - 1);
    1322          21 :   if (typ(VPOLGA) != t_VECSMALL) VPOLGA = gtovecsmall(VPOLGA);
    1323         112 :   for (i = 1; i < l; i++)
    1324             :   {
    1325          91 :     long e = VPOLGA[i];
    1326          91 :     if (e)
    1327             :     {
    1328          56 :       GEN D = divisorsu(i);
    1329          56 :       long j, lD = lg(D);
    1330         161 :       for (j = 1; j < lD; j++) v[ D[j] ] += e;
    1331             :     }
    1332             :   }
    1333         112 :   for (i = 1, cn = cd = 0; i < l; i++)
    1334             :   {
    1335          91 :     long e = v[i];
    1336          91 :     if (e > 0) cn += e; else cd -= e;
    1337             :   }
    1338          21 :   if (!cn || !cd) pari_err_TYPE("hgmgammatocyclo", VPOLGA);
    1339          14 :   *pE = n = cgetg(cn+1, t_VECSMALL);
    1340          14 :   *pD = d = cgetg(cd+1, t_VECSMALL);
    1341          98 :   for (i = cn = cd = 1; i < l; i++)
    1342             :   {
    1343          84 :     long e = v[i], j;
    1344         189 :     if (e < 0) for (j = 1; j <= -e; j++) d[cd++] = i;
    1345         126 :     else if (e > 0) for (j = 1; j <= e; j++) n[cn++] = i;
    1346             :   }
    1347          14 : }
    1348             : 
    1349             : static void
    1350          21 : hgmgammatoalpha(GEN VPOLGA, GEN *pA, GEN *pB)
    1351          21 : { hgmgammatocyclo(VPOLGA, pA, pB); hgmcyclotoalpha(pA, pB); }
    1352             : 
    1353             : /* A and B sorted */
    1354             : static long
    1355      554932 : zv_intersect(GEN A, GEN B)
    1356             : {
    1357      554932 :   long a = 1, b = 1, lA = lg(A), lB = lg(B);
    1358     1801464 :   while (a < lA && b < lB)
    1359             :   {
    1360     1697794 :     if      (A[a] < B[b]) a++;
    1361     1379231 :     else if (A[a] > B[b]) b++; else return 1;
    1362             :   }
    1363      103670 :   return 0;
    1364             : }
    1365             : static void
    1366         231 : remove_intersect(GEN *pA, GEN *pB)
    1367             : {
    1368         231 :   GEN V, W, A = *pA, B = *pB;
    1369         231 :   long a = 1, b = 1, v = 1, w = 1, lA, lB;
    1370         231 :   *pA = V = cgetg_copy(A, &lA);
    1371         231 :   *pB = W = cgetg_copy(B, &lB);
    1372        1267 :   while (a < lA && b < lB)
    1373             :   {
    1374        1036 :     int s = gcmp(gel(A,a), gel(B,b));
    1375        1036 :     if   (s < 0) gel(V,v++) = gel(A,a++);
    1376         707 :     else if (s > 0) gel(W,w++) = gel(B,b++);
    1377         273 :     else { a++; b++; }
    1378             :   }
    1379         364 :   while (a < lA) gel(V,v++) = gel(A,a++);
    1380         259 :   while (b < lB) gel(W,w++) = gel(B,b++);
    1381         231 :   setlg(V,v); setlg(W,w);
    1382         231 : }
    1383             : 
    1384             : /* al, be normalized, sorted, unknown intersection */
    1385             : static GEN
    1386         231 : albe2u(GEN al, GEN be)
    1387             : {
    1388             :   GEN ga;
    1389         231 :   remove_intersect(&al, &be);
    1390         231 :   ga = hgmalphatogamma(al, be);
    1391         231 :   return F2v_factorback(ga);
    1392             : }
    1393             : static GEN
    1394         469 : get_u(GEN al, GEN be, GEN CYCLOE, GEN VPOLGA, long DEG, long WT)
    1395             : {
    1396         469 :   GEN u, u0 = F2v_factorback(VPOLGA);
    1397         469 :   long b1 = get_b1(CYCLOE);
    1398         469 :   if (odd(DEG))
    1399             :   {
    1400         147 :     be = QV_normalize(RgV_addhalf(be));
    1401         147 :     u = albe2u(al, be);
    1402         147 :     if ((DEG & 3) == 3) u = negi(u);
    1403             :   }
    1404         322 :   else if (odd(WT)) { u = u0; if ((b1 & 3) == 2) u = negi(u); }
    1405             :   else
    1406             :   {
    1407          84 :     al = QV_normalize(RgV_addhalf(al));
    1408          84 :     u = shifti(albe2u(al, be), 1);
    1409          84 :     if (((DEG + b1) & 3) == 1) u = negi(u);
    1410             :   }
    1411         469 :   u0 = shifti(u0, 1); if ((b1 & 3) == 3) u0 = negi(u0);
    1412         469 :   return mkvec2(coredisc(u), u0);
    1413             : }
    1414             : 
    1415             : static long
    1416          35 : zv_sumeuler(GEN v)
    1417             : {
    1418          35 :   long i, l = lg(v);
    1419          35 :   GEN s = gen_0;
    1420         133 :   for (i = 1; i < l; i++)
    1421             :   {
    1422          98 :     if (v[i] <= 0) pari_err_TYPE("hgminit", v);
    1423          98 :     s = addiu(s, eulerphiu(v[i]));
    1424             :   }
    1425          35 :   return itou(s);
    1426             : }
    1427             : 
    1428             : /* (a, b): format (1) if denominator, format (2) if no denominator,
    1429             :  * format (3) if b not vector. */
    1430             : static GEN
    1431         518 : hgminit_i(GEN a, GEN b)
    1432             : {
    1433         518 :   long ta = typ(a), l = lg(a);
    1434         518 :   if (ta != t_VEC && ta != t_VECSMALL) pari_err_TYPE("hgminit", a);
    1435         511 :   if (!b)
    1436             :   {
    1437          84 :     if (l == 1) initab(a, a); /* error */
    1438          77 :     if (ta == t_VECSMALL || RgV_is_ZV(a))
    1439          49 :     {
    1440             :       long i;
    1441          56 :       if (ta != t_VECSMALL) a = vec_to_vecsmall(a);
    1442         161 :       for (i = 1; i < l; i++)
    1443         126 :         if (a[i] <= 0) break;
    1444          56 :       if (i != l)
    1445          21 :         hgmgammatoalpha(a, &a, &b); /* gamma */
    1446             :       else
    1447             :       { /* cyclo */
    1448          35 :         b = const_vecsmall(zv_sumeuler(a), 1);
    1449          35 :         hgmcyclotoalpha(&a, &b);
    1450             :       }
    1451             :     }
    1452             :     else /* alpha */
    1453          21 :       b = zerovec(l - 1);
    1454             :   }
    1455             :   else
    1456             :   {
    1457         427 :     if (typ(b) != ta) pari_err_TYPE("hgminit", b);
    1458         420 :     if (l > 1 && (ta == t_VECSMALL || (RgV_is_ZV(a) && RgV_is_ZV(b))))
    1459         259 :       hgmcyclotoalpha(&a, &b); /* cyclo */
    1460             :   }
    1461         483 :   return initab(a, b);
    1462             : }
    1463             : GEN
    1464         518 : hgminit(GEN val, GEN vbe)
    1465         518 : { pari_sp av = avma; return gerepilecopy(av, hgminit_i(val, vbe)); }
    1466             : 
    1467             : GEN
    1468          21 : hgmalpha(GEN hgm)
    1469             : {
    1470             :   GEN al, be;
    1471          21 :   if (!checkhgm(hgm)) pari_err_TYPE("hgmalpha", hgm);
    1472          14 :   al = hgm_get_VAL(hgm); be = hgm_get_VBE(hgm);
    1473          14 :   if (hgm_get_SWAP(hgm)) swap(al, be);
    1474          14 :   retmkvec2(gcopy(al), gcopy(be));
    1475             : }
    1476             : GEN
    1477          21 : hgmgamma(GEN hgm)
    1478             : {
    1479          21 :   pari_sp av = avma;
    1480             :   GEN g;
    1481          21 :   if (!checkhgm(hgm)) pari_err_TYPE("hgmgamma", hgm);
    1482          14 :   g = hgm_get_VPOLGA(hgm);
    1483          14 :   if (hgm_get_SWAP(hgm)) g = zv_neg(g);
    1484          14 :   return gerepilecopy(av, g);
    1485             : }
    1486             : GEN
    1487           7 : hgmcyclo(GEN hgm)
    1488             : {
    1489           7 :   pari_sp av = avma;
    1490             :   GEN A, B, C;
    1491           7 :   if (!checkhgm(hgm)) pari_err_TYPE("hgmcyclo", hgm);
    1492           7 :   C = hgm_get_CYCLOE(hgm); A = gel(C,1); B = gel(C,2);
    1493           7 :   if (hgm_get_SWAP(hgm)) swap(A, B);
    1494           7 :   return gerepilecopy(av, cycloE2cyclo(A, B));
    1495             : }
    1496             : 
    1497             : GEN
    1498          21 : hgmtwist(GEN hgm)
    1499             : {
    1500          21 :   pari_sp av = avma;
    1501             :   GEN val, vbe;
    1502          21 :   if (!checkhgm(hgm)) pari_err_TYPE("hgmtwist", hgm);
    1503          21 :   val = hgm_get_VAL(hgm);
    1504          21 :   vbe = hgm_get_VBE(hgm);
    1505          21 :   if (hgm_get_SWAP(hgm)) swap(val, vbe);
    1506          21 :   val = sort(RgV_addhalf(val));
    1507          21 :   vbe = sort(RgV_addhalf(vbe));
    1508          21 :   return gerepilecopy(av, initab(val, vbe));
    1509             : }
    1510             : 
    1511             : GEN
    1512         161 : hgmparams(GEN hgm)
    1513             : {
    1514         161 :   pari_sp av = avma;
    1515             :   GEN H, M;
    1516             :   long TT, DEG, WT;
    1517         161 :   if (!checkhgm(hgm)) pari_err_TYPE("hgmparams", hgm);
    1518         161 :   H = zx_to_ZX(hgm_get_HODGE(hgm));
    1519         161 :   TT = hgm_get_TT(hgm); DEG = hgm_get_DEG(hgm);
    1520         161 :   WT = hgm_get_WT(hgm); M = hgm_get_MVALUE(hgm);
    1521         161 :   return gerepilecopy(av, mkvec4(utoipos(DEG), utoi(WT),
    1522             :                                  mkvec2(H,stoi(TT)), M));
    1523             : }
    1524             : 
    1525             : /* symmetry at 1? */
    1526             : long
    1527           7 : hgmissymmetrical(GEN hgm)
    1528             : {
    1529             :   GEN A, B, C;
    1530             :   long a, i, j, lA, lB;
    1531           7 :   if (!checkhgm(hgm)) pari_err_TYPE("hgmissymmetrical", hgm);
    1532           7 :   C = hgm_get_CYCLOE(hgm);
    1533           7 :   if (odd(hgm_get_DEG(hgm))) return 0;
    1534           7 :   A = gel(C,1); lA = lg(A);
    1535           7 :   B = gel(C,2); lB = lg(B);
    1536          21 :   for (i = 1; i < lA; i++) if ((a = A[i]))
    1537             :   {
    1538           7 :     switch(i & 3)
    1539             :     { /* polcyclo(i)[-X] = polcyclo(j) */
    1540           0 :       case 0: j = i; break;
    1541           7 :       case 2: j = i >> 1; break;
    1542           0 :       default: j = i << 1; break;
    1543             :     }
    1544           7 :     if (j >= lB || B[j] != a) return 0;
    1545             :   }
    1546           7 :   return 1;
    1547             : }
    1548             : 
    1549             : /***************************************************************/
    1550             : /*         PART VI: Experimental euler factors                 */
    1551             : /***************************************************************/
    1552             : /* FIXME: one prime at a time */
    1553             : static GEN
    1554      227731 : hgmmodif(GEN an, GEN PPol)
    1555             : {
    1556      227731 :   pari_sp av = avma;
    1557      227731 :   long i, L = lg(an), lP = lg(PPol);
    1558             : 
    1559      674401 :   for (i = 1; i < lP; i++)
    1560             :   {
    1561      446670 :     GEN E = gel(PPol, i), pol = gel(E, 2);
    1562      446670 :     long d = degpol(pol);
    1563      446670 :     if (d)
    1564             :     {
    1565      433629 :       GEN v = vec_ei(L, 1);
    1566      433629 :       long j, p = itos(gel(E, 1)), pj = p;
    1567     1384271 :       for (j = 1; j <= d && pj <= L; j++, pj *= p)
    1568      950642 :         gel(v, pj) = RgX_coeff(pol, j);
    1569      433629 :       an = dirdiv(an, v);
    1570             :     }
    1571             :   }
    1572      227731 :   return gerepileupto(av, an);
    1573             : }
    1574             : 
    1575             : /***************************************************************/
    1576             : /*             PART VII: Make tables of HGMs                   */
    1577             : /***************************************************************/
    1578             : 
    1579             : static int
    1580         378 : ok_part(GEN v, GEN W)
    1581             : {
    1582         378 :   long l = lg(v), j;
    1583        1617 :   for (j = 1; j < l; j++)
    1584        1435 :     if (!gel(W,v[j])) return 0;
    1585         182 :   return 1;
    1586             : }
    1587             : static GEN
    1588          21 : mkphi()
    1589             : {
    1590          21 :   GEN W = const_vec(20, NULL);
    1591          21 :   gel(W,1) = mkvecsmall2(1,2);
    1592          21 :   gel(W,2) = mkvecsmall3(3,4,6);
    1593          21 :   gel(W,4) = mkvecsmall4(5,8,10,12);
    1594          21 :   gel(W,6) = mkvecsmall4(7,9,14,18);
    1595          21 :   gel(W,8) = mkvecsmall5(15,16,20,24,30);
    1596          21 :   gel(W,10)= mkvecsmall2(11,22);
    1597          21 :   gel(W,12)= mkvecsmalln(6, 13L,21L,26L,28L,36L,42L);
    1598          21 :   gel(W,16)= mkvecsmalln(6, 17L,32L,34L,40L,48L,60L);
    1599          21 :   gel(W,18)= mkvecsmall4(19,27,38,54);
    1600          21 :   gel(W,20)= mkvecsmall5(25,33,44,50,66);
    1601          21 :   return W;
    1602             : }
    1603             : 
    1604             : static GEN
    1605         182 : mkal(GEN p, GEN W)
    1606             : {
    1607             :   GEN res, V;
    1608         182 :   long l = lg(p), lV, i, j;
    1609         182 :   V = cgetg(l, t_VECSMALL);
    1610         987 :   for (i = 1; i < l; i++) V[i] = lg(gel(W, p[i])) - 1;
    1611         182 :   V = cyc2elts(V); lV = lg(V);
    1612         182 :   res = cgetg(1, t_VEC);
    1613       20769 :   for (j = 1; j < lV; j++)
    1614             :   {
    1615       20587 :     GEN t = cgetg(l, t_VECSMALL), v = gel(V,j);
    1616      157185 :     for (i = 1; i < l; i++) t[i] = umael2(W, p[i], v[i]+1);
    1617       20587 :     vecsmall_sort(t); if (!RgV_isin(res, t)) res = vec_append(res, t);
    1618             :   }
    1619         182 :   return res;
    1620             : }
    1621             : static GEN
    1622          21 : mkallal(long n)
    1623             : {
    1624          21 :   GEN W = mkphi(), p = partitions(n, NULL, NULL);
    1625          21 :   long i, c, l = lg(p);
    1626         399 :   for (i = c = 1; i < l; i++)
    1627         378 :     if (ok_part(gel(p,i), W)) gel(p, c++) = mkal(gel(p,i), W);
    1628          21 :   setlg(p, c); return shallowconcat1(p);
    1629             : }
    1630             : 
    1631             : static GEN
    1632          21 : mkalbe(long n)
    1633             : {
    1634          21 :   GEN w, L = mkallal(n);
    1635          21 :   long i, j, c, l = lg(L);
    1636          21 :   w = cgetg(l * (l / 2) + 1, t_VEC);
    1637        3941 :   for (i = c = 1; i < l; i++)
    1638             :   {
    1639        3920 :     GEN A = gel(L, i);
    1640        3920 :     long a = A[lg(A)-1];
    1641      558852 :     for (j = i + 1; j < l; j++)
    1642             :     {
    1643      554932 :       GEN B = gel(L, j);
    1644      554932 :       if (!zv_intersect(A, B))
    1645      179081 :         gel(w,c++) = (A[1]==1 || (B[1]!=1 && a > B[lg(B)-1]))?
    1646      179081 :                      mkvec2(B,A): mkvec2(A,B);
    1647             :     }
    1648             :   }
    1649          21 :   setlg(w,c); return w;
    1650             : }
    1651             : 
    1652             : static long
    1653      103670 : cyclowt(GEN a, GEN b)
    1654             : {
    1655      103670 :   pari_sp av = avma;
    1656             :   long TT;
    1657      103670 :   hgmcyclotoalpha(&a, &b);
    1658      103670 :   return gc_long(av, degpol(hodge(a, b, &TT)));
    1659             : }
    1660             : 
    1661             : GEN
    1662          21 : hgmbydegree(long n)
    1663             : {
    1664          21 :   pari_sp av = avma;
    1665          21 :   GEN w = cgetg(n+1, t_VEC), c = const_vecsmall(n, 1), v = mkalbe(n);
    1666          21 :   long i, l = lg(v);
    1667         154 :   for (i = 1; i <= n; i++) gel(w,i) = cgetg(l,t_VEC);
    1668      103691 :   for (i = 1; i < l; i++)
    1669             :   {
    1670      103670 :     GEN z = gel(v,i);
    1671      103670 :     long k = cyclowt(gel(z,1), gel(z,2)) + 1;
    1672      103670 :     gmael(w, k, c[k]++) = z;
    1673             :   }
    1674         154 :   for (i = 1; i <= n; i++) setlg(gel(w,i), c[i]);
    1675          21 :   return gerepilecopy(av, w);
    1676             : }
    1677             : 
    1678             : /***************************************************************/
    1679             : /*             PART VIII: L-function data                      */
    1680             : /***************************************************************/
    1681             : /* BAD sorted t_VECSMALL */
    1682             : static GEN
    1683         217 : removebad(GEN v, GEN BAD)
    1684             : {
    1685         217 :   long i, c, l = lg(v);
    1686         217 :   GEN w = cgetg(l, t_VECSMALL);
    1687         434 :   for (i = c = 1; i < lg(v); i++)
    1688         217 :     if (!zv_search(BAD, v[i])) w[c++] = v[i];
    1689         217 :   setlg(w, c); return w;
    1690             : }
    1691             : static GEN
    1692         609 : primedivisors(GEN t)
    1693         609 : { GEN fa = absZ_factor(t); return ZV_to_zv(gel(fa,1)); }
    1694             : 
    1695             : static GEN
    1696         385 : gacfac(long p, long m, long c)
    1697             : {
    1698         385 :   long i, l = p + m + 1;
    1699         385 :   GEN v = cgetg(l, t_VECSMALL);
    1700         728 :   for (i = 1; i <= p; i++) v[i] = -c;
    1701         658 :   for (     ; i <  l; i++) v[i] = 1 - c;
    1702         385 :   return v;
    1703             : }
    1704             : 
    1705             : static GEN
    1706         203 : hgmfindvga(GEN hgm, GEN t)
    1707             : {
    1708         203 :   GEN v, HODGE = hgm_get_HODGE(hgm);
    1709         203 :   long WT = degpol(HODGE), WT2 = (WT - 1) >> 1, i, c, h, fl = gequal1(t);
    1710         203 :   v = cgetg(WT + 2, t_VEC);
    1711         413 :   for (i = 0, c = 1; i <= WT2; i++)
    1712             :   {
    1713         210 :     if ((h = HODGE[i+2]))
    1714             :     {
    1715         210 :       if (fl && 2*i == WT - 1) h--;
    1716         210 :       if (h) gel(v, c++) = gacfac(h, h, i);
    1717             :     }
    1718             :   }
    1719         203 :   if (!odd(WT))
    1720             :   {
    1721         182 :     long h = HODGE[WT2+3], hp = h >> 1, hm;
    1722         182 :     if (!fl)
    1723             :     {
    1724         182 :       hm = hp;
    1725         182 :       if (odd(h))
    1726         154 :       { if (gsigne(t) >= 0 && gcmpgs(t, 1) <= 0) hm++; else hp++; }
    1727             :       else
    1728          28 :       { if (gcmpgs(t, 1) > 0) { hp++; hm--; } }
    1729         182 :       if (odd(hgm_get_TT(hgm) + WT2 + 1)) lswap(hp, hm);
    1730             :     }
    1731             :     else
    1732             :     {
    1733           0 :       h--; hm = h - hp;
    1734           0 :       if (odd(h) && odd(hgm_get_TT(hgm) + WT2 + 1)) lswap(hp, hm);
    1735             :     }
    1736         182 :     if (hp || hm) gel(v, c++) = gacfac(hp, hm, WT2 + 1);
    1737             :   }
    1738         203 :   if (c == 1) return cgetg(1, t_VECSMALL);
    1739         203 :   setlg(v, c); v = shallowconcat1(v);
    1740         203 :   vecsmall_sort(v); return v;
    1741             : }
    1742             : 
    1743             : /* Return [VGA, k, BAD, COND] where VGA as in lfun, k non motivic
    1744             :  * weight (s -> k - s), BAD is the list of wild primes and Euler factors,
    1745             :  * COND is the tame part of the conductor */
    1746             : static GEN
    1747         203 : hgmlfuninfty(GEN hgm, GEN t)
    1748             : {
    1749         203 :   pari_sp av = avma;
    1750         203 :   GEN VGA = hgmfindvga(hgm, t), T0, T1, v;
    1751         203 :   GEN COND, t1 = gsubgs(t, 1), BAD = primedivisors(hgm_get_BAD(hgm));
    1752         203 :   long k = hgm_get_WT(hgm) + 1, i, j, l;
    1753             : 
    1754         203 :   if (isintzero(t1)) T1 = cgetg(1, t_VECSMALL);
    1755             :   else
    1756             :   {
    1757         196 :     GEN fa = absZ_factor(numer_i(t1)), P = gel(fa,1), E = gel(fa,2);
    1758         196 :     if (!odd(k)) T1 = removebad(ZV_to_zv(P), BAD);
    1759             :     else
    1760             :     {
    1761         182 :       long j, lP = lg(P);
    1762         182 :       T1 = cgetg(lP, t_VECSMALL);
    1763         350 :       for (i = j = 1; i < lP; i++)
    1764             :       {
    1765         168 :         long p = itos(gel(P,i));
    1766         168 :         if (mpodd(gel(E,i)) && !zv_search(BAD, p)) T1[j++] = p;
    1767             :       }
    1768         182 :       setlg(T1,j);
    1769             :     }
    1770             :   }
    1771         203 :   COND = zv_prod_Z(T1);
    1772         203 :   T0 = removebad(gconcat(primedivisors(numer_i(t)),
    1773             :                          primedivisors(denom_i(t))), BAD);
    1774         217 :   for (i = 1; i < lg(T0); i++)
    1775             :   {
    1776          14 :     long p = T0[i], e = tameexpo(hgm, Q_lval(t, p));
    1777          14 :     COND = mulii(powuu(p, e), COND);
    1778             :   }
    1779         203 :   l = lg(BAD); v = cgetg(l, t_VEC);
    1780         490 :   for (i = j = 1; i < l; i++) /* FIXME: precompute wild Euler factors */
    1781         287 :     if (hgmclass(hgm, BAD[i], t) == C_BAD)
    1782         252 :       gel(v,j++) = mkvec2(utoipos(BAD[i]), pol_1(0));
    1783         203 :   setlg(v, j); return gerepilecopy(av, mkvec4(VGA, utoi(k), v, COND));
    1784             : }
    1785             : 
    1786             : /***************************************************************/
    1787             : /*             PART IX: GUESS OF WILD PRIMES                   */
    1788             : /***************************************************************/
    1789             : static GEN
    1790        1708 : vecmellin(GEN L, GEN K, GEN t, GEN td, GEN vroots, long bitprec)
    1791             : {
    1792        1708 :   long i, N = lfunthetacost(L, t, 0, bitprec);
    1793        1708 :   GEN v = cgetg(N+1, t_VEC);
    1794      427231 :   for (i = 1; i <= N; i++)
    1795      425523 :     gel(v,i) = gammamellininvrt(K, gmul(td, gel(vroots,i)), bitprec);
    1796        1708 :   return v;
    1797             : }
    1798             : 
    1799             : /* List of Weil polynomials for prime p of degree d and weight w */
    1800             : static GEN
    1801        1610 : listweil_i(ulong d, ulong p, ulong w)
    1802             : {
    1803             :   GEN V;
    1804        1610 :   if (d == 0) return mkvec(pol_1(0));
    1805        1225 :   if (odd(d))
    1806             :   {
    1807             :     GEN P;
    1808         644 :     if (odd(w)) return cgetg(1, t_VEC);
    1809         420 :     V = listweil_i(d - 1, p, w);
    1810         420 :     P = monomial(powuu(p, w / 2), 1, 0);
    1811         420 :     return shallowconcat(gmul(gsubsg(1, P), V), gmul(gaddsg(1, P), V));
    1812             :   }
    1813         581 :   if (d == 2)
    1814             :   {
    1815         553 :     long q = upowuu(p, w), s = usqrt(4 * q), i;
    1816         553 :     GEN Q = utoi(q);
    1817         553 :     V = cgetg(2*s + 3, t_VEC);
    1818        4844 :     for (i = 1; i <= 2*s + 1; i++) gel(V,i) = mkpoln(3, Q, stoi(1+s-i), gen_1);
    1819         553 :     gel(V, 2*s + 2) = mkpoln(3, negi(Q), gen_0, gen_1);
    1820         553 :     return V;
    1821             :   }
    1822          28 :   if (d == 4)
    1823             :   {
    1824          28 :     long q = upowuu(p, w), s = usqrt(16 * q), a, c, is2 = usqrt(4 * q);
    1825          28 :     double s2 = 2 * sqrt((double)q);
    1826          28 :     GEN Q2 = sqru(q), W, A, mA, tmp, Q;
    1827          28 :     V = cgetg(s + 3, t_VEC);
    1828         252 :     for (a = 0; a <= s; a++)
    1829             :     {
    1830         224 :       long b, m = ceil(a * s2) - 2 * q, M = ((a * a) >> 2) + 2 * q;
    1831         224 :       GEN AQ = stoi(a * q), mAQ = stoi(-a * q);
    1832         224 :       A = stoi(a); mA = stoi(-a);
    1833         224 :       W = cgetg(2*(M - m) + 3, t_VEC);
    1834        1876 :       for (b = m, c = 1; b <= M; b++)
    1835             :       {
    1836        1652 :         if (a) gel(W, c++) = mkpoln(5, Q2, mAQ, stoi(b), mA, gen_1);
    1837        1652 :         gel(W, c++) = mkpoln(5, Q2, AQ, stoi(b), A, gen_1);
    1838             :       }
    1839         224 :       setlg(W, c); gel(V, a + 1) = W;
    1840             :     }
    1841          28 :     W = cgetg(2 * is2 + 2, t_VEC);
    1842          28 :     tmp = mkpoln(3, stoi(-q), gen_0, gen_1);
    1843          28 :     Q = utoipos(q);
    1844         147 :     for (a = 0, c = 1; a <= is2; a++)
    1845             :     {
    1846         119 :       A = stoi(a); mA = stoi(-a);
    1847         119 :       if (a) gel(W, c++) = gmul(tmp, mkpoln(3, Q, mA, gen_1));
    1848         119 :       gel(W, c++) = gmul(tmp, mkpoln(3, Q, A, gen_1));
    1849             :     }
    1850          28 :     setlg(W, c); gel(V, s + 2) = W;
    1851          28 :     return shallowconcat1(V);
    1852             :   }
    1853           0 :   pari_err_IMPL("d > 5 in listweil");
    1854             :   return NULL; /* LCOV_EXCL_LINE */
    1855             : }
    1856             : 
    1857             : /* Same, weight <= w, by decreasing weight */
    1858             : static GEN
    1859         490 : listweilallw_i(ulong d, ulong p, ulong w)
    1860             : {
    1861         490 :   long i = d? w: 0;
    1862         490 :   GEN V = cgetg(i+2, t_VEC);
    1863        1680 :   for (; i >= 0; i--) gel(V,i+1) = listweil_i(d, p, i);
    1864         490 :   return shallowconcat1(V);
    1865             : }
    1866             : 
    1867             : static long
    1868      127344 : sumdeg(GEN W, GEN M)
    1869             : {
    1870      127344 :   long i, l = lg(M), s = 0;
    1871      375851 :   for (i = 1; i < l; i++) s += degpol(gmael(W,i,M[i]+1));
    1872      127344 :   return s;
    1873             : }
    1874             : 
    1875             : static GEN
    1876          35 : mul(GEN a, GEN b, ulong lim)
    1877             : {
    1878          35 :   long na = lg(a)-1, nb = lg(b)-1, i, j, n;
    1879          35 :   GEN c = cgetg(na * nb + 1, t_VECSMALL);
    1880         378 :   for (i = n = 1; i <= na; i++)
    1881        2597 :     for (j = 1; j <= nb; j++)
    1882             :     {
    1883        2254 :       ulong m = umuluu_or_0(a[i], b[j]);
    1884        2254 :       if (m && m <= lim) c[n++] = m;
    1885             :     }
    1886          35 :   setlg(c, n); return c;
    1887             : }
    1888             : static GEN
    1889          98 : listcond(GEN BAD, GEN achi, ulong min, ulong max)
    1890             : {
    1891          98 :   long i, j, l = lg(BAD);
    1892          98 :   GEN P = cgetg(l, t_VEC), z;
    1893         231 :   for (i = 1; i < l; i++)
    1894             :   {
    1895         133 :     long p = BAD[i], v = achi[i];
    1896         133 :     gel(P,i) = upowers_u(p, ulogint(max, p) - v, upowuu(p, v));
    1897             :   }
    1898         133 :   z = gel(P,1); for (i = 2; i < l; i++) z = mul(z, gel(P,i), max);
    1899          98 :   if (min > 1)
    1900             :   {
    1901           0 :     l = lg(z);
    1902           0 :     for (i = j = 1; i < l; i++)
    1903           0 :       if ((ulong)z[i] >= min) z[j++] = z[i];
    1904           0 :     setlg(z, j);
    1905             :   }
    1906          98 :   vecsmall_sort(z); return z;
    1907             : }
    1908             : 
    1909             : /* Artificial lfundiv by zeta(s - (k - 1) / 2). */
    1910             : static GEN
    1911          35 : lfundivraw(GEN L)
    1912             : {
    1913          35 :   long k = itos(ldata_get_k(L)), i;
    1914          35 :   GEN v = ldata_get_gammavec(L);
    1915          35 :   i = RgV_isin(v, stoi(-(k - 1) >> 2));
    1916          35 :   if (!i) pari_err_BUG("lfundivraw");
    1917          35 :   L = shallowcopy(L);
    1918          35 :   gel(L, 3) = vecsplice(v, i);
    1919          35 :   setlg(L, 7); return L;
    1920             : }
    1921             : 
    1922             : /* Compute forvec vectors from v, sorted by increasing total degree */
    1923             : static GEN
    1924          98 : forvecsort(GEN vF, GEN v)
    1925             : {
    1926          98 :   GEN w, E = cyc2elts(v);
    1927          98 :   long i, l = lg(E);
    1928          98 :   w = cgetg(l, t_VECSMALL);
    1929      127442 :   for (i = 1; i < l; i++) w[i] = sumdeg(vF, gel(E, i));
    1930          98 :   return vecpermute(E, vecsmall_indexsort(w)); /* by increasing total degree */
    1931             : }
    1932             : 
    1933             : static GEN
    1934          98 : lfunhgmwild(GEN L, GEN H, GEN t, GEN BAD, long pole, GEN hint, long bitprec)
    1935             : {
    1936             :   GEN v, K, t0, t0r, t0ir, t0i, t0k, N0, vM, vD, val, PPOL, vF, achi;
    1937          98 :   long d, lM, iN, iM, i, k, k2, prec = nbits2prec(bitprec), lB = lg(BAD);
    1938             :   long BADprod, limdeg;
    1939          98 :   ulong minN = 1, maxN = 2048;
    1940             : 
    1941          98 :   v = cgetg(lB, t_VECSMALL); PPOL = cgetg(lB, t_VEC);
    1942         231 :   for (i = 1; i < lB; i++)
    1943             :   {
    1944         133 :     v[i] = itou( gmael(BAD,i,1) );
    1945         133 :     gel(PPOL,i) = shallowcopy(gel(BAD,i));
    1946             :   }
    1947          98 :   BAD = v;
    1948          98 :   BADprod = zv_prod(BAD);
    1949          98 :   achi = get_achi(H, t, BAD);
    1950          98 :   if (pole) L = lfundivraw(L);
    1951          98 :   limdeg = d = ldata_get_degree(L);
    1952          98 :   N0 = ldata_get_conductor(L);
    1953          98 :   if (hint) switch(typ(hint))
    1954             :   {
    1955             :     long l;
    1956           0 :     case t_INT:
    1957           0 :       limdeg = itos(hint);
    1958           0 :       if (limdeg < 0 || limdeg > d) pari_err_TYPE("lfunhgm [bad hint]", hint);
    1959           0 :       break;
    1960           0 :     case t_VEC:
    1961           0 :       l = lg(hint);
    1962           0 :       if (l > 1 && l < 4)
    1963             :       {
    1964           0 :         GEN t = gel(hint,1), r;
    1965           0 :         if (typ(t) != t_INT || signe(t) <= 0)
    1966           0 :           pari_err_TYPE("lfunhgm [bad hint]", hint);
    1967           0 :         t = dvmdii(t, N0, &r);
    1968           0 :         if (r != gen_0)
    1969           0 :           pari_err_TYPE("lfunhgm [bad hint]", hint);
    1970           0 :         minN = maxN = itou(t);
    1971           0 :         if (l == 3)
    1972             :         {
    1973           0 :           t = gel(hint,2);
    1974           0 :           if (typ(t) != t_INT || signe(t) < 0 || cmpis(t, d) > 0)
    1975           0 :             pari_err_TYPE("lfunhgm [bad hint]", hint);
    1976           0 :           limdeg = itos(t);
    1977             :         }
    1978             :       }
    1979             :   }
    1980          98 :   k = itos(ldata_get_k(L)); k2 = (k-1) >> 1;
    1981          98 :   K = gammamellininvinit(ldata_get_gammavec(L), 0, bitprec + 32);
    1982          98 :   t0 = sstoQ(11, 10); t0i = ginv(t0); t0k = gpowgs(t0, k);
    1983          98 :   t0r = gpow(t0, sstoQ(2,d), prec); t0ir = ginv(t0r);
    1984             :   /* vF[i] list of Weil polynomials F for p = BAD[i],
    1985             :    * F = prod_j Fj, |reciprocal root of Fj|^2 = p^(w + 1 - nj)
    1986             :    * 2v_p(lc(F)) = deg F * (w + 1) - sum_j deg Fj * nj; */
    1987          98 :   vF = cgetg(lB, t_VEC);
    1988          98 :   vD = cgetg(lB, t_VEC); /* vD[k][l] = sum_j deg Fj * nj for F = vF[k][l] */
    1989          98 :   v = cgetg(lB, t_VECSMALL);
    1990         231 :   for (i = 1; i < lB; i++)
    1991             :   {
    1992         133 :     GEN W = cgetg(limdeg+2, t_VEC), D;
    1993         133 :     long p = BAD[i], j, lW;
    1994         623 :     for (j = 0; j <= limdeg; j++) gel(W,j+1) = listweilallw_i(j, p, k-1);
    1995         133 :     gel(vF, i) = W = shallowconcat1(W);
    1996         133 :     lW = lg(W); v[i] = lW-1;
    1997         133 :     gel(vD,i) = D = cgetg(lW, t_VEC);
    1998       10136 :     for (j = 1; j < lW; j++)
    1999             :     {
    2000       10003 :       GEN F = gel(W,j);
    2001       10003 :       D[j] = degpol(F) * k - 2 * Z_lval(leading_coeff(F), p);
    2002             :     }
    2003             :   }
    2004          98 :   vM = forvecsort(vF, v); lM = lg(vM);
    2005          98 :   if (DEBUGLEVEL) { err_printf(" lM = %ld ", lM); err_flush(); }
    2006          98 :   L = shallowcopy(L);
    2007          98 :   val = cgetg(lB, t_VECSMALL);
    2008             :   for(;;)
    2009           0 :   {
    2010             :     GEN z, vroots, an0, vN;
    2011             :     long lN, lim;
    2012             : 
    2013          98 :     z = listcond(BAD, achi, minN, maxN);
    2014          98 :     if (maxN == minN) /* from hint */
    2015             :     {
    2016           0 :       minN = 1; /* in case it fails */
    2017           0 :       maxN--;
    2018             :     }
    2019             :     else
    2020             :     {
    2021          98 :       minN = maxN+1;
    2022          98 :       maxN <<= 2;
    2023             :     }
    2024          98 :     lN = lg(z);
    2025          98 :     if (lN == 1) continue;
    2026          98 :     vN = cgetg(lN, t_VEC);
    2027        1785 :     for (i = 1; i < lN; i++) gel(vN, i) = muliu(N0, z[i]);
    2028          98 :     gel(L, 5) = gel(vN, lN - 1);
    2029          98 :     if (DEBUGLEVEL) err_printf("\nmaxN = %ld\n", itos(gel(L,5)));
    2030          98 :     lim = lfunthetacost(L, t0i, 0, bitprec);
    2031          98 :     vroots = dirpowers(lim, sstoQ(2, d), prec + EXTRAPREC64);
    2032          98 :     an0 = hgmcoefs(H, t, lim);
    2033          98 :     if (pole)
    2034             :     {
    2035          35 :       GEN w = vecpowuu(lim, k2);
    2036       44450 :       for (i = 1; i <= lim; i++)
    2037       44415 :         if (cgcd(i, BADprod) > 1) gel(w, i) = gen_0;
    2038          35 :       an0 = dirdiv(an0, w);
    2039             :     }
    2040         854 :     for (iN = 1; iN < lN; iN++)
    2041             :     {
    2042         854 :       pari_sp av0 = avma, av2;
    2043         854 :       GEN vecK, vecKi, N = gel(vN,iN), isqN = gpow(N, sstoQ(-1,d), prec);
    2044         854 :       if (DEBUGLEVEL) { err_printf(" %ld ", itos(N)); err_flush(); }
    2045         854 :       gel(L,5) = N;
    2046         854 :       vecK = vecmellin(L, K, t0, gmul(t0r,  isqN), vroots, bitprec);
    2047         854 :       vecKi= vecmellin(L, K, t0i,gmul(t0ir, isqN), vroots, bitprec);
    2048        2310 :       for (i = 1; i < lB; i++) val[i] = Z_lval(N, BAD[i]);
    2049         854 :       setlg(an0, lg(vecKi)); av2 = avma;
    2050     2305926 :       for (iM = 1; iM < lM; iM++, set_avma(av2))
    2051             :       {
    2052     2305170 :         GEN M = gel(vM, iM), an, eno, den;
    2053     3214589 :         for (i = 1; i < lB; i++)
    2054             :         {
    2055     2986858 :           GEN F = gmael(vF, i, M[i]+1);
    2056     2986858 :           long D, dF = degpol(F), a = val[i];
    2057     2986858 :           if (a < d - dF) break;
    2058     2797417 :           if (a)
    2059             :           {
    2060     2606261 :             if (d == dF) break;
    2061     1163771 :             D = mael(vD, i, M[i]+1); /* sum_j deg Fj nj */
    2062     1163771 :             if (D == d && d - dF != a) break;
    2063             :             /* a = a(chi) + d - deg F - 1 */
    2064     1100575 :             if (D == d-1 && a != d - dF + achi[i] - 1 && !gequal1(t)) break;
    2065             :           }
    2066      909419 :           gmael(PPOL, i, 2) = F;
    2067             :         }
    2068     2305170 :         if (i < lB) continue;
    2069      227731 :         an = hgmmodif(an0, PPOL);
    2070      227731 :         den = gmul(t0k, RgV_dotproduct(vecK, an));
    2071      227731 :         eno = gdiv(RgV_dotproduct(vecKi, an), den);
    2072      227731 :         if (gexpo(gsubgs(gabs(eno, prec), 1)) < -bitprec / 2)
    2073             :         {
    2074          98 :           if (pole)
    2075          77 :             for (i = 1; i < lB; i++)
    2076             :             {
    2077          42 :               GEN t = deg1pol_shallow(negi(powuu(BAD[i], k2)), gen_1, 0);
    2078          42 :               gmael(PPOL, i, 2) = gmul(gmael(PPOL, i, 2), t);
    2079             :             }
    2080         231 :           for (i = 1; i < lB; i++) gmael(PPOL, i, 2) = ginv(gmael(PPOL, i, 2));
    2081          98 :           eno = grndtoi(eno, NULL);
    2082          98 :           if (typ(eno) != t_INT) pari_err_BUG("lfunhgmwild");
    2083          98 :           return mkvec3(N, mkvec2(t, PPOL), eno);
    2084             :         }
    2085             :       }
    2086         756 :       set_avma(av0);
    2087             :     }
    2088             :   }
    2089             : }
    2090             : 
    2091             : static GEN
    2092          91 : BAD2small(GEN BAD)
    2093             : {
    2094          91 :   long i, l = lg(BAD);
    2095          91 :   GEN v = cgetg(l, t_VECSMALL);
    2096         210 :   for (i = 1; i < l; i++) v[i] = itou( gmael(BAD,i,1) );
    2097          91 :   return v;
    2098             : }
    2099             : 
    2100             : /* moments of motive */
    2101             : static GEN
    2102          91 : hgmmoments(GEN H, GEN t, GEN M, long nb)
    2103             : {
    2104          91 :   pari_sp av = avma, av2;
    2105          91 :   GEN P, v, S, L = hgmlfuninfty(H, t), BAD = BAD2small(gel(L, 3));
    2106          91 :   GEN k2 = gmul2n(gsubgs(gel(L,2), 1), -1);
    2107          91 :   long ct = 0, i, j, lP, lm, tm = typ(M);
    2108             : 
    2109          91 :   if (!nb) nb = 1000;
    2110          91 :   P = primes_zv(nb); lP = lg(P);
    2111          91 :   v = hgmcoefs(H, t, P[lP - 1]);
    2112          91 :   if (tm != t_VEC) M = gtovec(M);
    2113          91 :   av2 = avma; lm = lg(M);
    2114          91 :   S = const_vec(lm - 1, real_0(DEFAULTPREC));
    2115        3731 :   for (i = 1; i < lP; i++)
    2116             :   {
    2117        3640 :     long p = P[i];
    2118        3640 :     if (!Q_lval(t, p) && !zv_search(BAD, p))
    2119             :     {
    2120        3500 :       GEN T = gdiv(gel(v, p), gpow(utoipos(p), k2, DEFAULTPREC));
    2121        3500 :       ct++;
    2122        7000 :       for (j = 1; j < lm; j++)
    2123        3500 :         gel(S,j) = gadd(gel(S,j), gpow(T, gel(M,j), DEFAULTPREC));
    2124             :     }
    2125        3640 :     if ((i & 0xf) == 0) S = gerepilecopy(av2, S);
    2126             :   }
    2127          91 :   if (tm != t_VEC && tm != t_COL && tm != t_VECSMALL) S = gel(S, 1);
    2128          91 :   return gerepileupto(av, gdivgu(S, ct));
    2129             : }
    2130             : 
    2131             : /* Heuristic guess: is there a pole ? */
    2132             : static long
    2133          91 : lfunhgmispole(GEN H, GEN t, long nb)
    2134             : {
    2135          91 :   pari_sp av = avma;
    2136             :   GEN M;
    2137          91 :   if (!nb) nb = 40;
    2138          91 :   M = hgmmoments(H, t, gen_1, nb);
    2139          91 :   return gc_bool(av, gexpo(M) > -2);
    2140             : }
    2141             : 
    2142             : static GEN
    2143         112 : tag(GEN x, long t) { return mkvec2(mkvecsmall(t), x); }
    2144             : 
    2145             : static GEN
    2146         112 : lfunhgm_i(GEN hgm, GEN t, GEN hint, long bitprec)
    2147             : {
    2148         112 :   GEN L, vr, v = hgmlfuninfty(hgm, t), vga = zv_to_ZV(gel(v,1)), k = gel(v,2);
    2149         112 :   GEN BAD = gel(v,3), COND = gel(v, 4);
    2150         112 :   long pole = mpodd(k) && lfunhgmispole(hgm, t, 0), lB = lg(BAD);
    2151             : 
    2152         112 :   L = mkvecn(pole? 7: 6, tag(mkvec2(hgm,t), t_LFUN_HGM), gen_0, vga, k, COND,
    2153             :              gen_0, mkvec(mkvec2(shifti(addiu(k,1),-1), gen_0)));
    2154         112 :   if (pole && ldata_get_degree(L) == 1)
    2155             :   { /* motive = zeta */
    2156             :     long i;
    2157           0 :     gel(L,5) = gel(L,6) = gel(L,7) = gen_1;
    2158           0 :     if (lB != 1)
    2159             :     {
    2160           0 :       GEN R = mkrfrac(gen_1, deg1pol_shallow(gen_m1,gen_1,0));
    2161           0 :       gmael3(L, 1, 2, 2) = mkvec2(t, BAD);
    2162           0 :       for (i = 1; i < lB; i++) gmael(BAD, i, 2) = R;
    2163             :     }
    2164           0 :     return L;
    2165             :   }
    2166         112 :   if (lB == 1)
    2167             :   {
    2168          14 :     vr = lfunrootres(L, bitprec);
    2169          14 :     gel(L, 6) = gel(vr, 3);
    2170          14 :     if (pole) gel(L, 7) = gel(vr, 1);
    2171          14 :     return L;
    2172             :   }
    2173          98 :   v = lfunhgmwild(L, hgm, t, BAD, pole, hint, bitprec);
    2174          98 :   gel(L, 5) = gel(v, 1); /* N */
    2175          98 :   gmael3(L, 1, 2, 2) = gel(v, 2); /* [t, PPOL] */
    2176          98 :   gel(L, 6) = gel(v, 3); /* w */
    2177          98 :   if (pole)
    2178             :   {
    2179          35 :     vr = lfunrootres(L, bitprec);
    2180          35 :     gel(L, 7) = gel(vr, 1);
    2181             :   }
    2182          98 :   return L;
    2183             : }
    2184             : GEN
    2185         119 : lfunhgm(GEN hgm, GEN t, GEN hint, long bit)
    2186             : {
    2187         119 :   pari_sp av = avma;
    2188         119 :   if (!checkhgm(hgm)) pari_err_TYPE("lfunhgm", hgm);
    2189         112 :   return gerepilecopy(av, lfunhgm_i(hgm, t, hint, bit));
    2190             : }
    2191             : 
    2192             : GEN
    2193        3997 : dirhgm_worker(GEN P, ulong X, GEN hgm, GEN t)
    2194             : {
    2195        3997 :   pari_sp av = avma;
    2196        3997 :   long i, l = lg(P);
    2197        3997 :   GEN W = cgetg(l, t_VEC);
    2198       40821 :   for(i = 1; i < l; i++)
    2199             :   {
    2200       36824 :     ulong p = uel(P,i);
    2201       36824 :     long e, d = ulogint(X, p) + 1; /* minimal d such that p^d > X */
    2202       36823 :     gel(W,i) = RgXn_inv(hgmeulerfactorlimit(hgm, t, p, d-1, 0, &e), d);
    2203             :   }
    2204        3997 :   return gerepilecopy(av, mkvec2(P,W));
    2205             : }
    2206             : 
    2207             : GEN
    2208         399 : hgmcoefs(GEN hgm, GEN t, long n)
    2209             : {
    2210         399 :   GEN worker, bad = NULL;
    2211         399 :   if (!checkhgm(hgm)) pari_err_TYPE("hgmcoefs", hgm);
    2212         399 :   if (typ(t) == t_VEC && lg(t) == 3) { bad = gel(t,2); t = gel(t,1); }
    2213         399 :   if (!is_rational_t(typ(t))) pari_err_TYPE("hgmcoefs",t);
    2214         392 :   worker = snm_closure(is_entry("_dirhgm_worker"), mkvec2(hgm, t));
    2215         392 :   return pardireuler(worker, gen_2, stoi(n), NULL, bad);
    2216             : }

Generated by: LCOV version 1.16