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