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 : }
|
