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 318121 : pari_init_hgm(void) { HGMDATA2 = HGMDATA3 = NULL; }
25 :
26 : void
27 313185 : pari_close_hgm(void)
28 313185 : { 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 979398 : RgV_mirror(GEN x) { pari_APPLY_same(gsubsg(1, gel(x,i))); }
46 :
47 : static GEN
48 104139 : hodge(GEN val, GEN vbe, long *pTT)
49 : {
50 104139 : long r = lg(val) - 1, T, s, k;
51 104139 : GEN H, w = indexsort( shallowconcat(val, RgV_mirror(vbe)) );
52 1854657 : for (k = 1, s = T = 0; k <= 2*r; k++)
53 1750518 : if (w[k] > r) { if (--s < T) T = s; } else s++;
54 104139 : H = zero_zv(r + 1 - T);
55 1854657 : for (k = 1, s = 0; k <= 2*r; k++) /* (x^-T * Hodge) += x^s, s - T >= 0 */
56 1750518 : if (w[k] > r) s--; else { H[s + 1 - T]++; s++; }
57 104139 : *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 1400 : al2cyE(GEN v)
64 : {
65 : GEN dv, E, F, w;
66 1400 : long N, i, j, l = lg(v);
67 1400 : if (l == 1) return cgetg(1, t_VECSMALL);
68 1232 : w = Q_remove_denom(v, &dv);
69 1232 : if (!dv) return mkvecsmall(l-1);
70 987 : N = itou(dv); w = ZV_to_Flv(w, N); vecsmall_sort(w);
71 987 : E = zero_zv(N); /* to record polcyclo(i)^E[i] */
72 987 : F = cgetg(l, t_VECSMALL); /* to check input */
73 4879 : for (i = j = 1; i < l; i++)
74 : {
75 3892 : long k, m = w[i];
76 3892 : if (!m) { E[1]++; F[j++] = 0; }
77 3521 : else if (N % m == 0)
78 : {
79 1890 : long d = N/m, km;
80 1890 : E[d]++;
81 8946 : for (k = 1, km = m; k <= d; k++, km += m)
82 7056 : if (ugcd(k,d) == 1) F[j++] = km;
83 : }
84 : }
85 987 : setlg(F,j); vecsmall_sort(F);
86 987 : return gequal(w,F)? E: NULL;
87 : }
88 :
89 : static int
90 700 : cyE_intersect(GEN a, GEN b)
91 : {
92 700 : long i, l = minss(lg(a), lg(b));
93 2653 : for (i = 1; i < l; i++) if (a[i] && b[i]) return 1;
94 693 : return 0;
95 : }
96 : static GEN
97 700 : get_CYCLOE(GEN val, GEN vbe)
98 : {
99 700 : GEN Ea = al2cyE(val), Eb = al2cyE(vbe);
100 700 : if (!Ea || !Eb || cyE_intersect(Ea, Eb))
101 7 : pari_err_TYPE("hgminit [not a Q motive]", mkvec2(val,vbe));
102 693 : return mkvec2(Ea, Eb);
103 : }
104 :
105 : /* R[n/d] += mu(d) * r, for all d | n */
106 : static void
107 10472 : moebiusadd(GEN R, long n, long r)
108 : {
109 10472 : if (r) {
110 1764 : GEN D = divisorsu_moebius(gel(factoru(n), 1));
111 1764 : long j, l = lg(D);
112 1764 : R[n] += r; /* d = 1 */
113 3647 : for (j = 2; j < l; j++) { long d = D[j]; R[n / labs(d)] += d < 0? -r: r; }
114 : }
115 10472 : }
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 693 : get_VPOLGA(GEN E)
119 : {
120 693 : GEN Ea = gel(E,1), Eb = gel(E,2), R;
121 693 : long i, la = lg(Ea), lb = lg(Eb), n = maxuu(la, lb) - 1;
122 : pari_sp av;
123 693 : R = zero_zv(n); av = avma;
124 7728 : for (i = 1; i < la; i++) moebiusadd(R, i, Ea[i]);
125 4130 : for (i = 1; i < lb; i++) moebiusadd(R, i, -Eb[i]);
126 5726 : for (i = n; i; i--) if (R[i]) break;
127 693 : 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 1204 : QV_normalize(GEN a)
159 : {
160 1204 : long l = lg(a), i;
161 1204 : GEN v = cgetg(l, t_VEC);
162 5908 : for (i = 1; i < l; i++)
163 : {
164 4718 : GEN c = gel(a, i);
165 4718 : if (!is_rational_t(typ(c)))
166 14 : pari_err_TYPE("hgminit [not rational params]",c);
167 4704 : gel(v, i) = gfrac(c);
168 : }
169 1190 : return sort(v);
170 : }
171 :
172 : static GEN
173 1344 : mangoldtexp(long n)
174 : {
175 1344 : GEN D = divisorsu_moebius(gel(factoru(n),1)), P = gen_1;
176 1344 : long i, l = lg(D);
177 4081 : for (i = 1; i < l; i++)
178 : { /* prod (n/d) ^ (n/d mu(d)) */
179 2737 : long mud = D[i], nd = n / mud; /* mud = mu(d) * d */
180 2737 : P = gmul(P, powis(utoi(labs(nd)), nd));
181 : }
182 1344 : return P;
183 : }
184 : static GEN
185 924 : E2exp(GEN E)
186 : {
187 924 : long n, r, l = lg(E);
188 924 : GEN P = gen_1;
189 9786 : for (n = 1; n < l; n++)
190 8862 : if ((r = E[n])) P = gmul(P, gpowgs(mangoldtexp(n), r));
191 924 : return P;
192 : }
193 : static long
194 959 : 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 504 : initab(GEN a, GEN b)
200 : {
201 : GEN VALNUM, VALDEN, VBENUM, VBEDEN, CYCLOE, SIGNPAR;
202 : GEN BAD, HODGE, VPOLGA, res;
203 504 : long j, WT, TT, OFFMPOL, l = lg(a), DEG = l-1, SWAP = 0;
204 :
205 504 : if (lg(b) != l) pari_err_TYPE("hgminit [#al != #be]", mkvec2(a,b));
206 504 : if (l == 1) pari_err_TYPE("hgminit [#al = 0]", mkvec2(a,b));
207 490 : a = QV_normalize(a);
208 483 : b = QV_normalize(b);
209 476 : 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 469 : VALNUM = cgetg(l, t_VECSMALL);
216 469 : VALDEN = cgetg(l, t_VECSMALL);
217 469 : VBENUM = cgetg(l, t_VECSMALL);
218 469 : VBEDEN = cgetg(l, t_VECSMALL);
219 2443 : for (j = 1; j < l; j++)
220 : {
221 1974 : Qtoss(gel(a,j), &VALNUM[j], &VALDEN[j]);
222 1974 : Qtoss(gel(b,j), &VBENUM[j], &VBEDEN[j]);
223 : }
224 469 : BAD = negi(lcmii(ZV_lcm(zv_to_ZV(VALDEN)), ZV_lcm(zv_to_ZV(VBEDEN))));
225 469 : HODGE = hodge(a, b, &TT);
226 469 : WT = degpol(HODGE);
227 469 : CYCLOE = get_CYCLOE(a, b);
228 462 : VPOLGA = get_VPOLGA(CYCLOE);
229 462 : SIGNPAR = odd(WT)? gen_1: discprod(gel(CYCLOE, odd(DEG)? 2: 1));
230 462 : OFFMPOL = (get_b1(CYCLOE) - zv_sum(VPOLGA)) >> 1;
231 462 : res = cgetg(13, t_VEC);
232 462 : gel(res, 1) = a;
233 462 : gel(res, 2) = b;
234 462 : gel(res, 3) = CYCLOE;
235 462 : gel(res, 4) = VBENUM;
236 462 : gel(res, 5) = VBEDEN;
237 462 : gel(res, 6) = gdiv(E2exp(gel(CYCLOE,1)), E2exp(gel(CYCLOE,2))); /*MVALUE*/
238 462 : gel(res, 7) = VPOLGA;
239 462 : gel(res, 8) = BAD;
240 462 : gel(res, 9) = HODGE;
241 462 : gel(res, 10) = get_u(a, b, CYCLOE, VPOLGA, DEG, WT);
242 462 : gel(res, 11) = SIGNPAR;
243 462 : gel(res, 12) = mkvecsmall3(OFFMPOL, TT, SWAP);
244 462 : return res;
245 : }
246 : static int
247 1281 : checkhgm(GEN v)
248 : {
249 1260 : return typ(v) == t_VEC && lg(v) == 13
250 2541 : && typ(gel(v,12)) == t_VECSMALL && lg(gel(v,12)) == 4;
251 : }
252 :
253 : /* accessors */
254 141735 : 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 1582 : static GEN hgm_get_CYCLOE(GEN v) { return gel(v, 3); }
257 37596 : static GEN hgm_get_VBENUM(GEN v) { return gel(v, 4); }
258 37596 : static GEN hgm_get_VBEDEN(GEN v) { return gel(v, 5); }
259 623 : static GEN hgm_get_MVALUE(GEN v) { return gel(v, 6); }
260 687872 : static GEN hgm_get_VPOLGA(GEN v) { return gel(v, 7); }
261 37833 : static GEN hgm_get_BAD(GEN v) { return gel(v, 8); }
262 109716 : 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 141700 : static long hgm_get_DEG(GEN v) { return lg(hgm_get_VAL(v))-1; }
267 652492 : static long hgm_get_OFFMPOL(GEN v) { return gel(v, 12)[1]; }
268 76967 : static long hgm_get_TT(GEN v) { return gel(v, 12)[2]; }
269 672 : static long hgm_get_SWAP(GEN v) { return gel(v, 12)[3]; }
270 109353 : 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 59348 : for (i = 2; i < p; i++) gel(v, i+1) = gdivgu(gel(v, i), i);
285 2155 : gel(w, 1) = shallowcopy(v);
286 261110 : for (k = s = 1, t = p; k <= L; k++)
287 : {
288 258948 : gel(v, s) = gdivgu(gadd(gel(v, t), gel(v, s)), k + p - 1);
289 258954 : if (++s > p) { s = 1; gel(w, k/p + 1) = shallowcopy(v); }
290 258955 : if (++t > p) t = 1;
291 : }
292 2162 : 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 116949 : gel(v,1) = z = leafcopy(x);
302 116949 : setvalp(z, valp(z)+1);
303 625637 : for (i = 1; i < N; i++)
304 : {
305 508691 : gel(v, i+1) = z = gmul(z, gaddgs(x, i));
306 508654 : setvalp(z, valp(z) + 1);
307 : }
308 116946 : 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 742528 : for (k = 1; k <= N; k++) S = gadd(S, gmul(gel(C,k), gmael(v, k+1, m0+1)));
318 116945 : 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 203895 : if (x >= pD) x %= pD;
330 203895 : z = cgetg(5, t_PADIC);
331 203894 : z[1] = evalprecp(D) | evalvalp(v);
332 203893 : gel(z,2) = gp;
333 203893 : gel(z,3) = gpD;
334 203893 : 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 104482 : for (m = 1; m <= q; m += 2, set_avma(av2))
368 : {
369 102460 : if (f > 1) { m0 = m%p; m1 = m / p; M = m1 + m0 * pf1; } else M = m0 = m;
370 102460 : x = umultop(M, iQ, gp, gpD, D);
371 102459 : vres[m+1] = S = sumbinom(vall, x, m0, N, pD);
372 102458 : if (!odd(m0)) S = pD - S;
373 102458 : 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 10634 : for (v = 1, vpo = 2; vpo < Q; v++, vpo <<= 1)
387 : {
388 8613 : if (v == vva) continue;
389 93792 : for (m = vpo; m <= q; m += 2*vpo, set_avma(av2))
390 : {
391 86945 : ulong ms2 = (m >> 1) + 1;
392 86945 : 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 86945 : m0 = p - m%p; m1 = m/p + 1; M = pM + m1 - m0 * pf1;
396 86945 : mp = umultop(M, iQ, gp, gpD, D);
397 86945 : S = Rg_to_Fl(gmul(Qp_exp(gmul(mp, l2)), shifti(utoi(t), m0-1)), pD);
398 86943 : vres[m+1] = S;
399 86943 : if (odd(m0)) S = pD - S;
400 86943 : vres[pf - m] = Fl_inv(S, pD);
401 : }
402 : }
403 2021 : 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 89912 : 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 33936 : teichmodp2(GEN VPOLGA, ulong p, ulong p2)
475 : {
476 : ulong pN;
477 : long c, l, N;
478 33936 : GEN W = cgetg_copy(VPOLGA,&l); /* W[N] = teich(N)^N */
479 142872 : for (N = 1, pN = p; N < l; N++, pN += p)
480 108935 : if ((c = VPOLGA[N]))
481 : {
482 82251 : long N0; (void)z_lvalrem(N, p, &N0);
483 82251 : if (c < 0) N0 = Fl_inv(N0 % p, p);
484 82251 : W[N] = Fl_powu(N0, pN, p2);
485 : }
486 33937 : return W;
487 : }
488 : /* GF[i] = i! mod p^2 */
489 : static GEN
490 33934 : get_GF(ulong p, ulong p2)
491 : {
492 33934 : GEN F = cgetg(p, t_VECSMALL);
493 : ulong i;
494 33562885 : F[1] = 1; F[2] = 2; for (i = 3; i < p; i++) F[i] = Fl_mul(F[i-1], i, p2);
495 32181 : return F;
496 : }
497 :
498 : /* p odd, return v s.t. v[i] = 1/i mod p, i < p/2 */
499 : static GEN
500 33934 : Flv_inv_p2(ulong p)
501 : {
502 33934 : long i, g = pgener_Fl(p), ph = (p >> 1);
503 33936 : GEN v = cgetg(ph + 1, t_VECSMALL);
504 33936 : pari_sp av = avma;
505 33936 : GEN w = cgetg(ph, t_VECSMALL);
506 33367 : w[1] = g;
507 16862064 : for (i = 2; i < ph; i++) w[i] = Fl_mul(w[i-1], g, p); /* w[i]=g^i */
508 17375617 : for (i = 1; i < ph; i++)
509 : {
510 17341693 : long x = w[i], y = w[ph - i]; /* g^(-i) = - g^(ph - i) */
511 17341693 : if (x > ph) v[p - x] = y; else v[x] = p - y; /* 1/(-x) = y, 1/x = -y */
512 : }
513 33924 : 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 33933 : get_GH(GEN F, ulong p)
519 : {
520 33933 : long i, g, ph = p >> 1;
521 33933 : GEN H = cgetg(p, t_VECSMALL), v = Flv_inv_p2(p);
522 33936 : H[p-1] = (F[p-1] + 1) / p; g = p - H[p-1];
523 16310841 : for (i = 1; i <= ph; i++)
524 : {
525 16276918 : long c = g = Fl_add(g, v[i], p);
526 16278806 : H[i] = c = Fl_mul(c, i, p);
527 16280499 : H[p-1 - i] = Fl_neg(Fl_add(c, g, p), p);
528 : }
529 33923 : return H;
530 : }
531 : /* p odd, GPB[m+1] = Gamma_p(m/(p-1)) mod p^2 */
532 : static GEN
533 33934 : doprecompmodp2(ulong p, ulong p2)
534 : {
535 33934 : GEN F = get_GF(p, p2), H = get_GH(F, p), GPV = cgetg(p, t_VECSMALL);
536 : ulong r;
537 32806329 : for (r = 1; r < p; r++)
538 : {
539 32772401 : long t = Fl_mul(F[r], 1 + p * H[r], p2);
540 32764474 : GPV[p - r] = odd(r)? t: p2 - t;
541 : }
542 33928 : return GPV;
543 : }
544 :
545 : /***************************************************************/
546 : /* PART III: Motive initializations depending on p */
547 : /***************************************************************/
548 : static GEN
549 37596 : get_L1(GEN hgm, long PFM1, long f)
550 : {
551 37596 : GEN VBEDEN = hgm_get_VBEDEN(hgm);
552 37596 : GEN VBENUM = hgm_get_VBENUM(hgm), v;
553 37596 : long DEG = hgm_get_DEG(hgm), j;
554 37596 : v = const_vecsmall(PFM1, f * hgm_get_TT(hgm));
555 164024 : for (j = 1; j <= DEG; j++)
556 126427 : if (PFM1 % VBEDEN[j] == 0)
557 : {
558 118307 : long t = ((-(PFM1 / VBEDEN[j]) * VBENUM[j]) % PFM1) + 1;
559 118307 : if (t <= 0) t += PFM1;
560 118307 : v[t] -= f;
561 : }
562 37597 : 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 33963 : get_L0(GEN hgm, ulong PM1)
569 : {
570 33963 : GEN perm, vt, vc, VL0, VPOLGA = hgm_get_VPOLGA(hgm);
571 33962 : long w, S, c, d, j, N, m, l = lg(VPOLGA), D = (l * (l-1)) >> 1;
572 33962 : vt = cgetg(D+1, t_VECSMALL);
573 33962 : vc = cgetg(D+1, t_VECSMALL);
574 109086 : for (N = 2, c = 1; N < l; N++)
575 75125 : if (VPOLGA[N])
576 : {
577 : ulong Q;
578 : long q;
579 195481 : for (q = 1, Q = 0; q <= N; q++, Q += PM1, c++)
580 : {
581 147114 : vt[c] = ceildivuu(Q, N);
582 147113 : vc[c] = VPOLGA[N];
583 : }
584 : }
585 33961 : setlg(vt,c); setlg(vc,c); perm = vecsmall_indexsort(vt);
586 33965 : vt = vecsmallpermute(vt, perm);
587 33965 : vc = vecsmallpermute(vc, perm);
588 33965 : w = vt[1];
589 147137 : for (j = 2, d = 1; j < c; j++)
590 113172 : if (vt[j] == w) vc[d] += vc[j];
591 : else
592 : {
593 79289 : d++;
594 79289 : vt[d] = w = vt[j];
595 79289 : vc[d] = vc[j];
596 : }
597 33965 : d++; vt[d] = PM1; vc[d] = 0; /* sentinels */
598 33965 : VL0 = cgetg(PM1+1, t_VECSMALL); S = hgm_get_OFFMPOL(hgm);
599 147231 : for (j = 1; j < d; j++)
600 : {
601 34888700 : for (m = vt[j]; m < vt[j+1]; m++) VL0[m+1] = S;
602 113260 : S -= vc[j+1];
603 : }
604 33971 : return VL0;
605 : }
606 : /* 2) direct computation (single m), good when p << N log N */
607 : static long
608 614869 : L0(GEN hgm, ulong p, long f, long PFM1, long m)
609 : {
610 614869 : GEN VPOLGA = hgm_get_VPOLGA(hgm);
611 614869 : long e, S = hgm_get_OFFMPOL(hgm) * f, l = lg(VPOLGA);
612 : ulong pem;
613 1858442 : for (e = 0, pem = m; e < f; e++, pem *= p)
614 : {
615 1243573 : long N, s = 0;
616 : ulong Npem;
617 9691535 : for (N = 1, Npem = pem; N < l; N++, Npem += pem)
618 8447962 : if (VPOLGA[N]) s += VPOLGA[N] * (long)(Npem / PFM1);
619 1243573 : S -= s;
620 : }
621 614869 : 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 984444 : gapnpow(GEN T, long p, long f, long PFM1, long N, ulong Nmpe)
656 : {
657 : GEN Vr;
658 : long Nm, i, S;
659 :
660 984444 : 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 235082 : hgmC(GEN VPOLGA, GEN GPV, GEN TEICH, ulong p, long f, ulong PFM1, ulong m, long dfp)
673 : {
674 235082 : GEN r = gen_1;
675 235082 : long c, e, N, l = lg(VPOLGA);
676 : ulong Nmpe, mpe;
677 700175 : for (e = 0, mpe = m; e < f; e++, mpe = Fl_mul(mpe, p, PFM1))
678 3819376 : for (N = 1, Nmpe = mpe; N < l; N++, Nmpe = Fl_add(Nmpe, mpe, PFM1))
679 3354281 : if ((c = VPOLGA[N]))
680 : { /* Gamma_p(frac(N*m*p^e/(p^f-1))) */
681 1449533 : GEN z = utoi(GPV[Nmpe + 1]);
682 1449534 : 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 1449528 : r = gmul(r, gpowgs(z, c));
685 : }
686 235065 : if (typ(r) != t_PADIC) r = cvtop(r, utoipos(p), dfp);
687 235065 : return r;
688 : }
689 :
690 : /* Same modulo p^2, Wm[N] = teich(N^(sign(VPOLGA[N]) * m)) */
691 : static ulong
692 16421193 : hgmCmodp2(GEN VPOLGA, GEN Wm, GEN GPV, ulong PFM1, ulong m, ulong p2)
693 : {
694 16421193 : long l = lg(VPOLGA), c, N;
695 16421193 : ulong res = 1, Nm;
696 64881325 : for (N = 1, Nm = m; N < l; N++, Nm = Fl_add(Nm, m, PFM1))
697 48507223 : if ((c = VPOLGA[N]))
698 : {
699 : ulong z;
700 37099923 : if (c > 0) z = GPV[Nm + 1];
701 : else
702 : {
703 19064035 : c = -c;
704 19064035 : if (!Nm) z = 1;
705 : else
706 : {
707 19044630 : z = GPV[PFM1 + 1 - Nm];
708 19044630 : if (!odd(Nm)) z = p2 - z; /* GPV[Nm+1]^(-1) by reflection formula */
709 : }
710 : }
711 37740679 : if (N > 1) z = Fl_mul(z, Wm[N], p2);
712 37695138 : z = Fl_powu(z, c, p2);
713 : /* z = (GPV[Nm+1] * teich(N^m))^VPOLGA[N] */
714 37236641 : res = res == 1? z: Fl_mul(res, z, p2);
715 : }
716 16302575 : return res;
717 : }
718 :
719 : /***************************************************************/
720 : /* PART IV: Motive functions depending on p, t */
721 : /***************************************************************/
722 : static long
723 37837 : get_dfp(GEN hgm, long p, long f)
724 : {
725 37837 : long DEG = hgm_get_DEG(hgm), WT = hgm_get_WT(hgm);
726 37835 : long dfp = ceil(log((double)2*DEG) / log((double)p) + f * (WT * 0.5 ) + 1e-5);
727 37835 : if (p == 2) dfp++; /* FIXME: why ? */
728 37835 : 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 3660 : VL0 = (V && f == 1)? get_L0(hgm, PFM1): NULL;
741 3660 : 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 618607 : for (; i < lV; i++)
756 : {
757 614945 : long m = V? V[i]: i-1;
758 614945 : long L = VL0? VL0[m+1]: L0(hgm, p, f, PFM1, m);
759 614946 : long e = L + VL1[m+1];
760 614946 : if (!V && e >= dfp) c = gen_0;
761 : else
762 : { /* FIXME: dfp is fishy in Jordantame, don't trust if V = NULL */
763 235077 : pari_sp av = avma;
764 235077 : c = hgmC(VPOLGA, GPV, TEICH, p, f, PFM1, m, dfp);
765 235065 : setvalp(c, e); if (odd(L ^ l0)) c = gneg(c);
766 235066 : c = gerepileupto(av, padic_to_Q(c));
767 : }
768 614946 : gel(v, V? i: i+1) = c;
769 : }
770 3662 : return v;
771 : }
772 : /* Same mod p^2, f = 1 */
773 : static GEN
774 33935 : hgmCallmodp2(GEN hgm, ulong p)
775 : {
776 33935 : GEN C, GPV, VL0, VL1, W1, Wm, VPOLGA = hgm_get_VPOLGA(hgm);
777 33935 : long l = lg(VPOLGA), TT = hgm_get_TT(hgm);
778 33935 : ulong m, PFM1 = p - 1, p2 = p * p;
779 :
780 33935 : if (p & HIGHMASK) pari_err_OVERFLOW("hgmCallmodp2");
781 33928 : VL0 = get_L0(hgm, PFM1);
782 33936 : VL1 = get_L1(hgm, PFM1, 1);
783 33936 : W1 = teichmodp2(VPOLGA, p, p2);
784 33936 : Wm = shallowcopy(W1);
785 33934 : GPV = doprecompmodp2(p, p2);
786 33928 : C = cgetg(PFM1+2, t_VECSMALL); C[1] = evalvarn(0);
787 33931 : C[2] = TT > 1? 0: (TT == 1? p : 1); /* m = 0 */
788 32218580 : for (m = 1; m < PFM1; m++)
789 : {
790 32652912 : long e = VL0[m + 1] + VL1[m + 1], j;
791 : ulong c;
792 32652912 : if (e >= 2) c = 0;
793 : else
794 : {
795 16348019 : c = hgmCmodp2(VPOLGA, Wm, GPV, PFM1, m, p2);
796 16298059 : if (odd(VL0[m + 1] ^ VL0[1])) c = Fl_neg(c, p2);
797 15879243 : if (e == 1) c = (c % p) * p;
798 : }
799 32184136 : C[m + 2] = c;
800 95618303 : for (j = 2; j < l; j++)
801 63433654 : 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 3850 : inv(ulong q, long n)
809 : {
810 3850 : ulong z = q + 1;
811 : long i;
812 7980 : for (i = 1; i < n; i++) z = z * q + 1;
813 3850 : return z;
814 : }
815 :
816 : /* General H function: C(teich(f + O(p^dfp))) / (1 - p^f) */
817 : static GEN
818 3850 : hgmH(GEN C, long p, long f, long dfp, GEN t)
819 : {
820 3850 : GEN q = powuu(p, dfp), z;
821 : long n;
822 3850 : (void)Q_lvalrem(t, p, &t);
823 3850 : z = Rg_to_Fp(t, q);
824 3850 : z = Zp_teichmuller(z, utoipos(p), dfp, q);
825 3850 : z = FpX_eval(C, z, q);
826 3850 : n = dfp / f; if (!(dfp % f)) n--;
827 3850 : z = Fp_mulu(z, inv(upowuu(p,f), n), q);
828 3850 : return Fp_center(z, q, shifti(q,-1));
829 : }
830 : static GEN
831 33923 : hgmHmodp2(GEN C, ulong p, GEN t)
832 : {
833 33923 : ulong p2 = p * p, wt = Fl_powu(Rg_to_Fl(t, p2), p, p2);
834 33933 : 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 37795 : hgmtrace(GEN hgm, long p, long f, GEN t, long c)
842 : {
843 37795 : long dfp = get_dfp(hgm, p, f);
844 37792 : if (c == C_FAKE) return hgmU(hgm, p, f, t, dfp);
845 37561 : 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 693 : F2v_factorback(GEN E)
851 : {
852 693 : long i, l = lg(E);
853 693 : GEN c = gen_1;
854 4186 : for (i = 1; i < l; i++)
855 3493 : if (odd(E[i])) c = muliu(c, i);
856 693 : 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 33934 : hgmsign(GEN hgm, long p, GEN t)
863 : {
864 : GEN u;
865 33934 : if (odd(hgm_get_WT(hgm))) return 1;
866 28418 : t = ginv(t); u = gmul(gsubsg(1, t), hgm_get_SIGNPAR(hgm));
867 28418 : 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 245 : hgmcharcond(GEN hgm, GEN t)
872 : {
873 : GEN u;
874 245 : 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 364 : upowers_u(ulong r, long n, ulong a)
893 : {
894 364 : long i, l = n + 2, t = a;
895 364 : GEN v = cgetg(l, t_VECSMALL);
896 2310 : for (i = 1; i < l; i++) { v[i] = t; t *= r; }
897 364 : return v;
898 : }
899 : static GEN
900 36875 : powers_u(ulong r, long n, ulong a)
901 : {
902 36875 : long i, l = n + 2;
903 36875 : GEN v = cgetg(l, t_VEC), t = utoi(a);
904 106589 : for (i = 1; i < l; i++) { gel(v,i) = t; t = muliu(t, r); }
905 36874 : return v;
906 : }
907 : static GEN
908 36874 : mkpowers(long p, long d, long WT)
909 : {
910 : ulong q, r;
911 36874 : if (WT == 0) q = r = 1;
912 35698 : 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 36874 : 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 39709 : 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 39498 : for (; k >= nE; k--) gel(T, k) = gen_0;
932 108607 : for (; k > 0; k--) gel(T, k) = gel(E, k+1);
933 36874 : return RgV_to_RgX(T,0);
934 : }
935 :
936 : static long
937 37630 : hgmclass(GEN hgm, long p, GEN t)
938 : {
939 37630 : GEN BAD = hgm_get_BAD(hgm);
940 : long ap, bp;
941 37630 : if (!umodiu(BAD, p))
942 : {
943 665 : long v = Q_lval(t, p);
944 665 : if (v && v + Q_lval(hgm_get_MVALUE(hgm), p) == 0)
945 : {
946 147 : GEN N = hgmcharcond(hgm, t);
947 147 : if (umodiu(N, p)) return C_FAKE; /* a(chi) = 0 */
948 : }
949 553 : return C_BAD;
950 : }
951 36964 : if (typ(t) == t_INT)
952 : {
953 23993 : 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 36949 : frobpoltrunc(GEN hgm, GEN t, long c, long p, long B, long *pF)
968 : {
969 36949 : long DEG = hgm_get_DEG(hgm), WT = hgm_get_WT(hgm);
970 36949 : long D = isint1(t)? (odd(WT) ? DEG - 2 : DEG - 1): DEG;
971 36949 : long f, mi, m, q = upowuu(p, WT >> 1);
972 : GEN E, s;
973 :
974 36948 : mi = minss(B, (c == C_FAKE)? D: D >> 1);
975 36948 : m = (mi == D && c == C_TAME1 && !odd(DEG))? mi: mi+1;
976 36948 : s = cgetg(m+1, t_POL); s[1] = evalsigne(1)|evalvarn(0);
977 74730 : for (f = 1; f < m; f++) gel(s, f+1) = negi( hgmtrace(hgm, p, f, t, c) );
978 36936 : *pF = 0; s = RgXn_expint(s, m);
979 36937 : if (mi == D) return s;
980 36874 : if (c == C_TAME1)
981 : {
982 2940 : long SIGN = kroiu(hgm_get_U(hgm), p);
983 2940 : if (odd(DEG))
984 203 : E = Efuneq(s, p, DEG-1, WT, SIGN, B);
985 : else
986 : {
987 2737 : GEN T = deg1pol_shallow(stoi(- SIGN * q), gen_1, 0);
988 2737 : E = RgXn_mul(s, RgXn_inv(T, m), m);
989 2737 : E = Efuneq(E, p, DEG - 2, WT, 1, B);
990 2737 : if (!gequal1(t) || !odd(WT)) E = gmul(E, T);
991 : }
992 2940 : *pF = 1;
993 2940 : if (!odd(WT))
994 : {
995 : GEN T, u, t0;
996 217 : long v = Q_lvalrem(gsubgs(t, 1), p, &t0);
997 217 : if (!odd(v))
998 : {
999 63 : if (typ(t0) == t_FRAC) t0 = mulii(gel(t0,1), gel(t0,2));
1000 63 : u = coredisc(mulii(t0, hgm_get_U0(hgm)));
1001 63 : T = deg1pol_shallow(stoi(-kroiu(u, p) * q), gen_1, 0);
1002 63 : E = gmul(E, T); *pF = 0;
1003 : }
1004 : }
1005 : }
1006 : else
1007 : {
1008 33934 : long SIGN = hgmsign(hgm, p, t);
1009 33934 : E = Efuneq(s, p, DEG, WT, SIGN, B);
1010 : }
1011 36874 : return E;
1012 : }
1013 :
1014 : GEN
1015 84 : hgmcoef(GEN hgm, GEN t, GEN n)
1016 : {
1017 84 : pari_sp av = avma;
1018 84 : GEN T, P, E, F = check_arith_all(n, "hgmcoef");
1019 : long i, lP;
1020 :
1021 84 : if (!checkhgm(hgm)) pari_err_TYPE("hgmcoef", hgm);
1022 84 : if (!is_rational_t(typ(t))) pari_err_TYPE("hgmcoef",t);
1023 77 : if (hgm_get_SWAP(hgm)) t = ginv(t);
1024 77 : if (!F) { F = Z_factor(n); P = gel(F,1); }
1025 : else
1026 : {
1027 21 : P = gel(F,1);
1028 21 : if (lg(P) == 1 || signe(gel(P,1)) <= 0) return gen_0;
1029 21 : n = typ(n) == t_VEC? gel(n,1): factorback(F);
1030 : }
1031 77 : if (signe(n) <= 0) pari_err_DOMAIN("hgmcoef", "n", "<=", gen_0, n);
1032 77 : E = gel(F,2); lP = lg(P); T = gen_1;
1033 77 : T = gen_1;
1034 161 : for (i = 1; i < lP; i++)
1035 : {
1036 84 : long e, p = itos(gel(P, i)), f = itos(gel(E, i)), c = hgmclass(hgm, p, t);
1037 : GEN A;
1038 84 : if (c == C_BAD) pari_err_IMPL("hgmcoef for bad primes");
1039 84 : A = frobpoltrunc(hgm, t, c, p, f, &e);
1040 84 : T = gmul(T, RgX_coeff(RgXn_inv(A, f+1), f));
1041 : }
1042 77 : return gerepilecopy(av, T);
1043 : }
1044 :
1045 : static GEN
1046 14 : count2list(GEN E)
1047 : {
1048 14 : long i, j, r, l = lg(E);
1049 14 : GEN v = cgetg(zv_sum(E)+1, t_VECSMALL);
1050 49 : for (i = j = 1; i < l; i++)
1051 56 : if ((r = E[i])) while(r--) v[j++] = i;
1052 14 : setlg(v,j); return v;
1053 : }
1054 : static GEN hgminit_i(GEN val, GEN vbe);
1055 : #if 0
1056 : static long
1057 : minval(GEN F, long p)
1058 : {
1059 : long i, d = degpol(F), m = LONG_MAX;
1060 : for (i = 1; i <= d; i++) m = minss(m, Z_lval(gel(F,i+2), p) / i);
1061 : return m;
1062 : }
1063 : static GEN
1064 : cycloE_filter(GEN A, long p)
1065 : {
1066 : long i, j, l = lg(A);
1067 : GEN B = shallowcopy(A);
1068 : if (p < l)
1069 : for (i = p; i < l; i += p)
1070 : if (A[i])
1071 : {
1072 : (void)z_lvalrem(i / p, p, &j);
1073 : B[j] += A[i];
1074 : B[i] = 0;
1075 : }
1076 : return B;
1077 : }
1078 : /* doesn't work because A/B is not valid hypergeometric data */
1079 : static GEN
1080 : eulfacbadnew(GEN hgm, GEN t, long p, long *pe)
1081 : {
1082 : GEN AB = hgm_get_CYCLOE(hgm), gp = utoipos(p);
1083 : GEN A = cycloE_filter(gel(AB,1), p);
1084 : GEN B = cycloE_filter(gel(AB,2), p);
1085 : GEN hgmp = hgminit_i(count2list(A), count2list(B));
1086 : GEN cF, E, F = hgmeulerfactor(hgmp, t, p, &E);
1087 : long s, v, d = degpol(F), w = hgm_get_WT(hgm);
1088 : *pe = 0;
1089 : if (!d) return pol_1(0);
1090 : cF = leading_coeff(F);
1091 : v = Z_lvalrem(cF, p, &cF);
1092 : if (!is_pm1(cF)) return pol_1(0);
1093 : s = minval(F, p);
1094 : if (s) { F = RgX_unscale(F, powis(gp, -s)); v -= s * d; }
1095 : if ((2 * v) % d || (!odd(w) && v % d)) return F;
1096 : return RgX_unscale(F, powis(gp, w / 2 - v / d));
1097 : }
1098 : #endif
1099 : static GEN Jordantameexpo(GEN hgm, long v, GEN t0, long p, long *pe);
1100 : static GEN
1101 37259 : hgmeulerfactorlimit(GEN hgm, GEN t, long p, long d, long flag, long *pe)
1102 : {
1103 37259 : long c = hgmclass(hgm, p, t);
1104 37257 : if (c == C_TAME0)
1105 : {
1106 91 : long v = Q_lvalrem(t, p, &t);
1107 91 : return Jordantameexpo(hgm, v, t, p, pe);
1108 : }
1109 37166 : if (c != C_BAD) return frobpoltrunc(hgm, t, c, p, d, pe);
1110 301 : if (flag) { *pe = -1; return gen_0; } else { *pe = 0; return pol_1(0); }
1111 : }
1112 :
1113 : GEN
1114 441 : hgmeulerfactor(GEN hgm, GEN t, long p, GEN* pE)
1115 : {
1116 441 : pari_sp av = avma;
1117 : long e, B;
1118 : GEN P;
1119 441 : if (!checkhgm(hgm)) pari_err_TYPE("hgmeulerfactor", hgm);
1120 441 : if (!is_rational_t(typ(t))) pari_err_TYPE("hgmeulerfactor",t);
1121 434 : if (hgm_get_SWAP(hgm)) t = ginv(t);
1122 434 : B = (long)(hgm_get_DEG(hgm) * log(p)) + 1;
1123 434 : P = gerepilecopy(av, hgmeulerfactorlimit(hgm, t, p, B, 1, &e));
1124 420 : if (pE) *pE = stoi(e);
1125 420 : return P;
1126 : }
1127 :
1128 : /***********************************************************************/
1129 : /* Tame Euler factors */
1130 : /***********************************************************************/
1131 : /* FIXME: implement properly like RgXn_sqrt */
1132 : static GEN
1133 42 : RgXn_sqrtnu(GEN h, long f, long e)
1134 : {
1135 42 : if (f == 1) return h;
1136 7 : if (f == 2) return RgXn_sqrt(h, e);
1137 0 : return ser2rfrac_i(gsqrtn(RgX_to_ser(h, e + 2), utoipos(f), NULL, 0));
1138 : }
1139 : static GEN
1140 77 : Jordantame(GEN hgm, GEN t0, long m, long p)
1141 : {
1142 : GEN P, T, C, ZP, V;
1143 : long d, phim, f, j, c, q, qm, dfp;
1144 :
1145 77 : if (m == 1)
1146 : {
1147 35 : long e = hgm_get_WT(hgm) - get_b1(hgm_get_CYCLOE(hgm));
1148 35 : return deg1pol_shallow(negi(powuu(p, (e+1) >> 1)), gen_1, 0);
1149 : }
1150 42 : phim = eulerphiu(m); f = Fl_order(p % m, phim, m);
1151 42 : d = phim + 1; /* DEG >= phim >= f */
1152 42 : q = upowuu(p, f); qm = (q - 1) / m;
1153 42 : V = cgetg(phim + 1, t_VECSMALL);
1154 133 : for (j = c = 1; j < m; j++)
1155 91 : if (cgcd(j, m) == 1) V[c++] = j * qm;
1156 42 : dfp = get_dfp(hgm, p, f);
1157 42 : C = hgmCall(hgm, p, f, dfp, V);
1158 42 : ZP = zeropadic_shallow(utoipos(p), dfp);
1159 42 : T = teich(gadd(t0, ZP)); P = pol_1(0);
1160 133 : for (j = 1; j < lg(V); j++)
1161 : {
1162 91 : GEN Q, c = gmul(gpowgs(T, V[j]), gel(C,j));
1163 91 : Q = RgXn_red_shallow(RgX_shift_shallow(RgX_Rg_mul(P, c), f), d);
1164 91 : P = RgX_sub(P, Q); /* P * (1 - c x^f) mod X^d */
1165 : }
1166 42 : return centerlift(RgXn_sqrtnu(P, f, d));
1167 : }
1168 :
1169 : static GEN
1170 105 : eulfactameinit(GEN hgm, long v)
1171 : {
1172 105 : GEN C = hgm_get_CYCLOE(hgm);
1173 105 : if (!v) pari_err_BUG("hgmeulerfactor [incorrect t in eulfactame]");
1174 105 : if (hgm_get_SWAP(hgm)) v = -v;
1175 105 : return gel(C, (v < 0)? 2: 1);
1176 : }
1177 : static long
1178 14 : tameexpo(GEN hgm, long v)
1179 : {
1180 14 : GEN W = eulfactameinit(hgm,v);
1181 14 : long e, m, l = lg(W);
1182 42 : for (m = 1, e = 0; m < l; m++)
1183 28 : if (W[m] && v % m == 0) e += eulerphiu(m);
1184 14 : return hgm_get_DEG(hgm) - e;
1185 : }
1186 : static GEN
1187 91 : Jordantameexpo(GEN hgm, long v, GEN t0, long p, long *pe)
1188 : {
1189 91 : GEN P = pol_1(0), W = eulfactameinit(hgm, v);
1190 91 : long e, m, l = lg(W);
1191 308 : for (m = 1, e = 0; m < l; m++)
1192 217 : if (W[m] && v % m == 0)
1193 : {
1194 77 : P = gmul(P, Jordantame(hgm, t0, m, p));
1195 77 : e += eulerphiu(m);
1196 : }
1197 91 : *pe = hgm_get_DEG(hgm) - e; return P;
1198 : }
1199 :
1200 : /***************************************************************/
1201 : /* PART IV.5: Fake wild primes for t */
1202 : /***************************************************************/
1203 : /* Compute g_q(r)=pi^sq(r)*gauss(r) */
1204 : static GEN
1205 3612 : gkgauss(long p, long f, GEN vp, long r, GEN ZP, GEN *sp)
1206 : {
1207 3612 : GEN S = gen_0, P = gen_m1;
1208 3612 : long i, qm1 = vp[f+1] - 1;
1209 13342 : for (i = 1; i <= f; i++)
1210 : {
1211 9730 : GEN t = gfrac(sstoQ(r * vp[i], qm1));
1212 9730 : S = gadd(S, t);
1213 9730 : P = gmul(P, Qp_gamma(gadd(t, ZP)));
1214 : }
1215 3612 : *sp = gmulsg(p - 1, S); /* t_INT */
1216 3612 : return P;
1217 : }
1218 :
1219 : static GEN
1220 1428 : hgmG(GEN hgm, long p, long f, GEN vp, long r, GEN ZP)
1221 : {
1222 1428 : GEN S = gen_0, P = gen_1, VPOLGA = hgm_get_VPOLGA(hgm);
1223 : long n, c;
1224 5040 : for (n = 1; n < lg(VPOLGA); n++)
1225 3612 : if ((c = VPOLGA[n]))
1226 : {
1227 3612 : GEN sq, g = gkgauss(p, f, vp, r*n, ZP, &sq);
1228 3612 : S = addii(S, mulsi(c, sq));
1229 3612 : P = gmul(P, gpowgs(g, c));
1230 : }
1231 : /* p - 1 | S */
1232 1428 : return gmul(P, powis(utoineg(p), itos(S) / (p - 1)));
1233 : }
1234 :
1235 : /* multiplicity of -r / (q-1) in beta */
1236 : static long
1237 2856 : hgmmulti(GEN B, long q, long r)
1238 : {
1239 2856 : long d = (q-1) / ugcd(r, q-1);
1240 2856 : return d >= lg(B)? 0: B[d];
1241 : }
1242 :
1243 : /* We must have w(M)^r*hgmQ(q,r)=hgmC(q,r)/hgmC(q,0) */
1244 : static GEN
1245 1428 : hgmQ(GEN hgm, long p, long f, GEN vp, long r, GEN ZP)
1246 : {
1247 1428 : pari_sp av = avma;
1248 1428 : GEN B = gel(hgm_get_CYCLOE(hgm), 2);
1249 1428 : long q = vp[f+1], m0 = hgmmulti(B, q, 0), m1 = hgmmulti(B, q, r);
1250 1428 : GEN c = powis(utoipos(q), hgm_get_TT(hgm) + m0 - m1);
1251 1428 : if (odd(m0)) c = negi(c);
1252 1428 : return gerepileupto(av, padic_to_Q(gmul(c, hgmG(hgm, p, f, vp, r, ZP))));
1253 : }
1254 :
1255 : static GEN
1256 231 : hgmU(GEN hgm, long p, long f, GEN t, long dfp)
1257 : {
1258 231 : pari_sp av = avma;
1259 231 : GEN ZP = zeropadic_shallow(utoipos(p), dfp), vp = upowers_u(p, f, 1);
1260 231 : long q = vp[f+1], i;
1261 231 : GEN Q = cgetg(q+1, t_POL);
1262 231 : Q[1] = evalsigne(1)|evalvarn(0);
1263 1659 : for (i = 2; i <= q; i++) gel(Q, i) = hgmQ(hgm, p, f, vp, i, ZP);
1264 231 : t = p == 2? gen_1: gmul(hgm_get_MVALUE(hgm), t);
1265 231 : return gerepileupto(av, hgmH(Q, p, f, dfp, t));
1266 : }
1267 :
1268 : /***************************************************************/
1269 : /* PART V: Utility programs and main init */
1270 : /***************************************************************/
1271 : static GEN
1272 7 : cycloE2cyclo(GEN A, GEN B)
1273 7 : { return mkvec2(count2list(A), count2list(B)); }
1274 : #if 0
1275 : GEN
1276 : hgmalphatocyclo(GEN val, GEN vbe)
1277 : { GEN C = get_CYCLOE(val,vbe); return cycloE2cyclo(gel(C,1), gel(C,2)); }
1278 : #endif
1279 :
1280 : static GEN
1281 843850 : allprims(long n, GEN cache)
1282 : {
1283 : long l, i, c;
1284 : GEN z, v;
1285 843850 : if (gel(cache,n)) return gel(cache,n);
1286 418635 : z = coprimes_zv(n); l = lg(z); v = cgetg(l, t_VEC);
1287 3036782 : for (i = c = 1; i < l; i++)
1288 2618147 : if (z[i]) gel(v, c++) = mkfracss(i, n);
1289 418635 : setlg(v, c); return gel(cache,n) = v;
1290 : }
1291 : static GEN
1292 207956 : zv_to_prims(GEN D, GEN cache)
1293 : {
1294 207956 : long i, l = lg(D);
1295 207956 : GEN v = cgetg(l, t_VEC);
1296 1051806 : for (i = 1; i < l; i++)
1297 : {
1298 843857 : if (D[i] <= 0) pari_err_TYPE("hgmcyclotoalpha", D);
1299 843850 : gel(v, i) = allprims(D[i], cache);
1300 : }
1301 207949 : return shallowconcat1(v);
1302 : }
1303 : static void
1304 103978 : hgmcyclotoalpha(GEN *pA, GEN *pB)
1305 : {
1306 103978 : GEN v, D = *pA, E = *pB;
1307 103978 : if (typ(D) != t_VECSMALL) D = gtovecsmall(D);
1308 103978 : if (typ(E) != t_VECSMALL) E = gtovecsmall(E);
1309 103978 : v = const_vec(maxss(vecsmall_max(D), vecsmall_max(E)), NULL);
1310 103978 : gel(v,1) = mkvec(gen_0);
1311 103978 : *pA = zv_to_prims(D, v);
1312 103978 : *pB = zv_to_prims(E, v);
1313 103971 : if (lg(*pA) != lg(*pB))
1314 0 : pari_err_TYPE("hgminit [incorrect lengths]", mkvec2(D,E));
1315 103971 : }
1316 :
1317 : static GEN
1318 231 : hgmalphatogamma(GEN val, GEN vbe) { return get_VPOLGA(get_CYCLOE(val, vbe)); }
1319 :
1320 : static void
1321 21 : hgmgammatocyclo(GEN VPOLGA, GEN *pD, GEN *pE)
1322 : {
1323 21 : long i, cn, cd, l = lg(VPOLGA);
1324 21 : GEN d, n, v = zero_zv(l - 1);
1325 21 : if (typ(VPOLGA) != t_VECSMALL) VPOLGA = gtovecsmall(VPOLGA);
1326 112 : for (i = 1; i < l; i++)
1327 : {
1328 91 : long e = VPOLGA[i];
1329 91 : if (e)
1330 : {
1331 56 : GEN D = divisorsu(i);
1332 56 : long j, lD = lg(D);
1333 161 : for (j = 1; j < lD; j++) v[ D[j] ] += e;
1334 : }
1335 : }
1336 112 : for (i = 1, cn = cd = 0; i < l; i++)
1337 : {
1338 91 : long e = v[i];
1339 91 : if (e > 0) cn += e; else cd -= e;
1340 : }
1341 21 : if (!cn || !cd) pari_err_TYPE("hgmgammatocyclo", VPOLGA);
1342 14 : *pE = n = cgetg(cn+1, t_VECSMALL);
1343 14 : *pD = d = cgetg(cd+1, t_VECSMALL);
1344 98 : for (i = cn = cd = 1; i < l; i++)
1345 : {
1346 84 : long e = v[i], j;
1347 189 : if (e < 0) for (j = 1; j <= -e; j++) d[cd++] = i;
1348 126 : else if (e > 0) for (j = 1; j <= e; j++) n[cn++] = i;
1349 : }
1350 14 : }
1351 :
1352 : static void
1353 21 : hgmgammatoalpha(GEN VPOLGA, GEN *pA, GEN *pB)
1354 21 : { hgmgammatocyclo(VPOLGA, pA, pB); hgmcyclotoalpha(pA, pB); }
1355 :
1356 : /* A and B sorted */
1357 : static long
1358 554932 : zv_intersect(GEN A, GEN B)
1359 : {
1360 554932 : long a = 1, b = 1, lA = lg(A), lB = lg(B);
1361 1801464 : while (a < lA && b < lB)
1362 : {
1363 1697794 : if (A[a] < B[b]) a++;
1364 1379231 : else if (A[a] > B[b]) b++; else return 1;
1365 : }
1366 103670 : return 0;
1367 : }
1368 : static void
1369 231 : remove_intersect(GEN *pA, GEN *pB)
1370 : {
1371 231 : GEN V, W, A = *pA, B = *pB;
1372 231 : long a = 1, b = 1, v = 1, w = 1, lA, lB;
1373 231 : *pA = V = cgetg_copy(A, &lA);
1374 231 : *pB = W = cgetg_copy(B, &lB);
1375 1267 : while (a < lA && b < lB)
1376 : {
1377 1036 : int s = gcmp(gel(A,a), gel(B,b));
1378 1036 : if (s < 0) gel(V,v++) = gel(A,a++);
1379 707 : else if (s > 0) gel(W,w++) = gel(B,b++);
1380 273 : else { a++; b++; }
1381 : }
1382 364 : while (a < lA) gel(V,v++) = gel(A,a++);
1383 259 : while (b < lB) gel(W,w++) = gel(B,b++);
1384 231 : setlg(V,v); setlg(W,w);
1385 231 : }
1386 :
1387 : /* al, be normalized, sorted, unknown intersection */
1388 : static GEN
1389 231 : albe2u(GEN al, GEN be)
1390 : {
1391 : GEN ga;
1392 231 : remove_intersect(&al, &be);
1393 231 : ga = hgmalphatogamma(al, be);
1394 231 : return F2v_factorback(ga);
1395 : }
1396 : static GEN
1397 462 : get_u(GEN al, GEN be, GEN CYCLOE, GEN VPOLGA, long DEG, long WT)
1398 : {
1399 462 : GEN u, u0 = F2v_factorback(VPOLGA);
1400 462 : long b1 = get_b1(CYCLOE);
1401 462 : if (odd(DEG))
1402 : {
1403 147 : be = QV_normalize(RgV_addhalf(be));
1404 147 : u = albe2u(al, be);
1405 147 : if ((DEG & 3) == 3) u = negi(u);
1406 : }
1407 315 : else if (odd(WT)) { u = u0; if ((b1 & 3) == 2) u = negi(u); }
1408 : else
1409 : {
1410 84 : al = QV_normalize(RgV_addhalf(al));
1411 84 : u = shifti(albe2u(al, be), 1);
1412 84 : if (((DEG + b1) & 3) == 1) u = negi(u);
1413 : }
1414 462 : u0 = shifti(u0, 1); if ((b1 & 3) == 3) u0 = negi(u0);
1415 462 : return mkvec2(coredisc(u), u0);
1416 : }
1417 :
1418 : static long
1419 35 : zv_sumeuler(GEN v)
1420 : {
1421 35 : long i, l = lg(v);
1422 35 : GEN s = gen_0;
1423 133 : for (i = 1; i < l; i++)
1424 : {
1425 98 : if (v[i] <= 0) pari_err_TYPE("hgminit", v);
1426 98 : s = addiu(s, eulerphiu(v[i]));
1427 : }
1428 35 : return itou(s);
1429 : }
1430 :
1431 : /* (a, b): format (1) if denominator, format (2) if no denominator,
1432 : * format (3) if b not vector. */
1433 : static GEN
1434 511 : hgminit_i(GEN a, GEN b)
1435 : {
1436 511 : long ta = typ(a), l = lg(a);
1437 511 : if (ta != t_VEC && ta != t_VECSMALL) pari_err_TYPE("hgminit", a);
1438 504 : if (!b)
1439 : {
1440 84 : if (l == 1) initab(a, a); /* error */
1441 77 : if (ta == t_VECSMALL || RgV_is_ZV(a))
1442 49 : {
1443 : long i;
1444 56 : if (ta != t_VECSMALL) a = vec_to_vecsmall(a);
1445 161 : for (i = 1; i < l; i++)
1446 126 : if (a[i] <= 0) break;
1447 56 : if (i != l)
1448 21 : hgmgammatoalpha(a, &a, &b); /* gamma */
1449 : else
1450 : { /* cyclo */
1451 35 : b = const_vecsmall(zv_sumeuler(a), 1);
1452 35 : hgmcyclotoalpha(&a, &b);
1453 : }
1454 : }
1455 : else /* alpha */
1456 21 : b = zerovec(l - 1);
1457 : }
1458 : else
1459 : {
1460 420 : if (typ(b) != ta) pari_err_TYPE("hgminit", b);
1461 413 : if (l > 1 && (ta == t_VECSMALL || (RgV_is_ZV(a) && RgV_is_ZV(b))))
1462 259 : hgmcyclotoalpha(&a, &b); /* cyclo */
1463 : }
1464 476 : return initab(a, b);
1465 : }
1466 : GEN
1467 511 : hgminit(GEN val, GEN vbe)
1468 511 : { pari_sp av = avma; return gerepilecopy(av, hgminit_i(val, vbe)); }
1469 :
1470 : GEN
1471 21 : hgmalpha(GEN hgm)
1472 : {
1473 : GEN al, be;
1474 21 : if (!checkhgm(hgm)) pari_err_TYPE("hgmalpha", hgm);
1475 14 : al = hgm_get_VAL(hgm); be = hgm_get_VBE(hgm);
1476 14 : if (hgm_get_SWAP(hgm)) swap(al, be);
1477 14 : retmkvec2(gcopy(al), gcopy(be));
1478 : }
1479 : GEN
1480 21 : hgmgamma(GEN hgm)
1481 : {
1482 21 : pari_sp av = avma;
1483 : GEN g;
1484 21 : if (!checkhgm(hgm)) pari_err_TYPE("hgmgamma", hgm);
1485 14 : g = hgm_get_VPOLGA(hgm);
1486 14 : if (hgm_get_SWAP(hgm)) g = zv_neg(g);
1487 14 : return gerepilecopy(av, g);
1488 : }
1489 : GEN
1490 7 : hgmcyclo(GEN hgm)
1491 : {
1492 7 : pari_sp av = avma;
1493 : GEN A, B, C;
1494 7 : if (!checkhgm(hgm)) pari_err_TYPE("hgmcyclo", hgm);
1495 7 : C = hgm_get_CYCLOE(hgm); A = gel(C,1); B = gel(C,2);
1496 7 : if (hgm_get_SWAP(hgm)) swap(A, B);
1497 7 : return gerepilecopy(av, cycloE2cyclo(A, B));
1498 : }
1499 :
1500 : GEN
1501 21 : hgmtwist(GEN hgm)
1502 : {
1503 21 : pari_sp av = avma;
1504 : GEN val, vbe;
1505 21 : if (!checkhgm(hgm)) pari_err_TYPE("hgmtwist", hgm);
1506 21 : val = hgm_get_VAL(hgm);
1507 21 : vbe = hgm_get_VBE(hgm);
1508 21 : if (hgm_get_SWAP(hgm)) swap(val, vbe);
1509 21 : val = sort(RgV_addhalf(val));
1510 21 : vbe = sort(RgV_addhalf(vbe));
1511 21 : return gerepilecopy(av, initab(val, vbe));
1512 : }
1513 :
1514 : GEN
1515 161 : hgmparams(GEN hgm)
1516 : {
1517 161 : pari_sp av = avma;
1518 : GEN H, M;
1519 : long TT, DEG, WT;
1520 161 : if (!checkhgm(hgm)) pari_err_TYPE("hgmparams", hgm);
1521 161 : H = zx_to_ZX(hgm_get_HODGE(hgm));
1522 161 : TT = hgm_get_TT(hgm); DEG = hgm_get_DEG(hgm);
1523 161 : WT = hgm_get_WT(hgm); M = hgm_get_MVALUE(hgm);
1524 161 : return gerepilecopy(av, mkvec4(utoipos(DEG), utoi(WT),
1525 : mkvec2(H,stoi(TT)), M));
1526 : }
1527 :
1528 : /* symmetry at 1? */
1529 : long
1530 7 : hgmissymmetrical(GEN hgm)
1531 : {
1532 : GEN A, B, C;
1533 : long a, i, j, lA, lB;
1534 7 : if (!checkhgm(hgm)) pari_err_TYPE("hgmissymmetrical", hgm);
1535 7 : C = hgm_get_CYCLOE(hgm);
1536 7 : if (odd(hgm_get_DEG(hgm))) return 0;
1537 7 : A = gel(C,1); lA = lg(A);
1538 7 : B = gel(C,2); lB = lg(B);
1539 21 : for (i = 1; i < lA; i++) if ((a = A[i]))
1540 : {
1541 7 : switch(i & 3)
1542 : { /* polcyclo(i)[-X] = polcyclo(j) */
1543 0 : case 0: j = i; break;
1544 7 : case 2: j = i >> 1; break;
1545 0 : default: j = i << 1; break;
1546 : }
1547 7 : if (j >= lB || B[j] != a) return 0;
1548 : }
1549 7 : return 1;
1550 : }
1551 :
1552 : /***************************************************************/
1553 : /* PART VI: Experimental euler factors */
1554 : /***************************************************************/
1555 : /* FIXME: one prime at a time */
1556 : static GEN
1557 227731 : hgmmodif(GEN an, GEN PPol)
1558 : {
1559 227731 : pari_sp av = avma;
1560 227731 : long i, L = lg(an), lP = lg(PPol);
1561 :
1562 674401 : for (i = 1; i < lP; i++)
1563 : {
1564 446670 : GEN E = gel(PPol, i), pol = gel(E, 2);
1565 446670 : long d = degpol(pol);
1566 446670 : if (d)
1567 : {
1568 433629 : GEN v = vec_ei(L, 1);
1569 433629 : long j, p = itos(gel(E, 1)), pj = p;
1570 1384271 : for (j = 1; j <= d && pj <= L; j++, pj *= p)
1571 950642 : gel(v, pj) = RgX_coeff(pol, j);
1572 433629 : an = dirdiv(an, v);
1573 : }
1574 : }
1575 227731 : return gerepileupto(av, an);
1576 : }
1577 :
1578 : /***************************************************************/
1579 : /* PART VII: Make tables of HGMs */
1580 : /***************************************************************/
1581 :
1582 : static int
1583 378 : ok_part(GEN v, GEN W)
1584 : {
1585 378 : long l = lg(v), j;
1586 1617 : for (j = 1; j < l; j++)
1587 1435 : if (!gel(W,v[j])) return 0;
1588 182 : return 1;
1589 : }
1590 : static GEN
1591 21 : mkphi()
1592 : {
1593 21 : GEN W = const_vec(20, NULL);
1594 21 : gel(W,1) = mkvecsmall2(1,2);
1595 21 : gel(W,2) = mkvecsmall3(3,4,6);
1596 21 : gel(W,4) = mkvecsmall4(5,8,10,12);
1597 21 : gel(W,6) = mkvecsmall4(7,9,14,18);
1598 21 : gel(W,8) = mkvecsmall5(15,16,20,24,30);
1599 21 : gel(W,10)= mkvecsmall2(11,22);
1600 21 : gel(W,12)= mkvecsmalln(6, 13L,21L,26L,28L,36L,42L);
1601 21 : gel(W,16)= mkvecsmalln(6, 17L,32L,34L,40L,48L,60L);
1602 21 : gel(W,18)= mkvecsmall4(19,27,38,54);
1603 21 : gel(W,20)= mkvecsmall5(25,33,44,50,66);
1604 21 : return W;
1605 : }
1606 :
1607 : static GEN
1608 182 : mkal(GEN p, GEN W)
1609 : {
1610 : GEN res, V;
1611 182 : long l = lg(p), lV, i, j;
1612 182 : V = cgetg(l, t_VECSMALL);
1613 987 : for (i = 1; i < l; i++) V[i] = lg(gel(W, p[i])) - 1;
1614 182 : V = cyc2elts(V); lV = lg(V);
1615 182 : res = cgetg(1, t_VEC);
1616 20769 : for (j = 1; j < lV; j++)
1617 : {
1618 20587 : GEN t = cgetg(l, t_VECSMALL), v = gel(V,j);
1619 157185 : for (i = 1; i < l; i++) t[i] = umael2(W, p[i], v[i]+1);
1620 20587 : vecsmall_sort(t); if (!RgV_isin(res, t)) res = vec_append(res, t);
1621 : }
1622 182 : return res;
1623 : }
1624 : static GEN
1625 21 : mkallal(long n)
1626 : {
1627 21 : GEN W = mkphi(), p = partitions(n, NULL, NULL);
1628 21 : long i, c, l = lg(p);
1629 399 : for (i = c = 1; i < l; i++)
1630 378 : if (ok_part(gel(p,i), W)) gel(p, c++) = mkal(gel(p,i), W);
1631 21 : setlg(p, c); return shallowconcat1(p);
1632 : }
1633 :
1634 : static GEN
1635 21 : mkalbe(long n)
1636 : {
1637 21 : GEN w, L = mkallal(n);
1638 21 : long i, j, c, l = lg(L);
1639 21 : w = cgetg(l * (l / 2) + 1, t_VEC);
1640 3941 : for (i = c = 1; i < l; i++)
1641 : {
1642 3920 : GEN A = gel(L, i);
1643 3920 : long a = A[lg(A)-1];
1644 558852 : for (j = i + 1; j < l; j++)
1645 : {
1646 554932 : GEN B = gel(L, j);
1647 554932 : if (!zv_intersect(A, B))
1648 179081 : gel(w,c++) = (A[1]==1 || (B[1]!=1 && a > B[lg(B)-1]))?
1649 179081 : mkvec2(B,A): mkvec2(A,B);
1650 : }
1651 : }
1652 21 : setlg(w,c); return w;
1653 : }
1654 :
1655 : static long
1656 103670 : cyclowt(GEN a, GEN b)
1657 : {
1658 103670 : pari_sp av = avma;
1659 : long TT;
1660 103670 : hgmcyclotoalpha(&a, &b);
1661 103670 : return gc_long(av, degpol(hodge(a, b, &TT)));
1662 : }
1663 :
1664 : GEN
1665 21 : hgmbydegree(long n)
1666 : {
1667 21 : pari_sp av = avma;
1668 21 : GEN w = cgetg(n+1, t_VEC), c = const_vecsmall(n, 1), v = mkalbe(n);
1669 21 : long i, l = lg(v);
1670 154 : for (i = 1; i <= n; i++) gel(w,i) = cgetg(l,t_VEC);
1671 103691 : for (i = 1; i < l; i++)
1672 : {
1673 103670 : GEN z = gel(v,i);
1674 103670 : long k = cyclowt(gel(z,1), gel(z,2)) + 1;
1675 103670 : gmael(w, k, c[k]++) = z;
1676 : }
1677 154 : for (i = 1; i <= n; i++) setlg(gel(w,i), c[i]);
1678 21 : return gerepilecopy(av, w);
1679 : }
1680 :
1681 : /***************************************************************/
1682 : /* PART VIII: L-function data */
1683 : /***************************************************************/
1684 : /* BAD sorted t_VECSMALL */
1685 : static GEN
1686 217 : removebad(GEN v, GEN BAD)
1687 : {
1688 217 : long i, c, l = lg(v);
1689 217 : GEN w = cgetg(l, t_VECSMALL);
1690 434 : for (i = c = 1; i < lg(v); i++)
1691 217 : if (!zv_search(BAD, v[i])) w[c++] = v[i];
1692 217 : setlg(w, c); return w;
1693 : }
1694 : static GEN
1695 609 : primedivisors(GEN t)
1696 609 : { GEN fa = absZ_factor(t); return ZV_to_zv(gel(fa,1)); }
1697 :
1698 : static GEN
1699 385 : gacfac(long p, long m, long c)
1700 : {
1701 385 : long i, l = p + m + 1;
1702 385 : GEN v = cgetg(l, t_VECSMALL);
1703 728 : for (i = 1; i <= p; i++) v[i] = -c;
1704 658 : for ( ; i < l; i++) v[i] = 1 - c;
1705 385 : return v;
1706 : }
1707 :
1708 : static GEN
1709 203 : hgmfindvga(GEN hgm, GEN t)
1710 : {
1711 203 : GEN v, HODGE = hgm_get_HODGE(hgm);
1712 203 : long WT = degpol(HODGE), WT2 = (WT - 1) >> 1, i, c, h, fl = gequal1(t);
1713 203 : v = cgetg(WT + 2, t_VEC);
1714 413 : for (i = 0, c = 1; i <= WT2; i++)
1715 : {
1716 210 : if ((h = HODGE[i+2]))
1717 : {
1718 210 : if (fl && 2*i == WT - 1) h--;
1719 210 : if (h) gel(v, c++) = gacfac(h, h, i);
1720 : }
1721 : }
1722 203 : if (!odd(WT))
1723 : {
1724 182 : long h = HODGE[WT2+3], hp = h >> 1, hm;
1725 182 : if (!fl)
1726 : {
1727 182 : hm = hp;
1728 182 : if (odd(h))
1729 154 : { if (gsigne(t) >= 0 && gcmpgs(t, 1) <= 0) hm++; else hp++; }
1730 : else
1731 28 : { if (gcmpgs(t, 1) > 0) { hp++; hm--; } }
1732 182 : if (odd(hgm_get_TT(hgm) + WT2 + 1)) lswap(hp, hm);
1733 : }
1734 : else
1735 : {
1736 0 : h--; hm = h - hp;
1737 0 : if (odd(h) && odd(hgm_get_TT(hgm) + WT2 + 1)) lswap(hp, hm);
1738 : }
1739 182 : if (hp || hm) gel(v, c++) = gacfac(hp, hm, WT2 + 1);
1740 : }
1741 203 : if (c == 1) return cgetg(1, t_VECSMALL);
1742 203 : setlg(v, c); v = shallowconcat1(v);
1743 203 : vecsmall_sort(v); return v;
1744 : }
1745 :
1746 : /* Return [VGA, k, BAD, COND] where VGA as in lfun, k non motivic
1747 : * weight (s -> k - s), BAD is the list of wild primes and Euler factors,
1748 : * COND is the tame part of the conductor */
1749 : static GEN
1750 203 : hgmlfuninfty(GEN hgm, GEN t)
1751 : {
1752 203 : pari_sp av = avma;
1753 203 : GEN VGA = hgmfindvga(hgm, t), T0, T1, v;
1754 203 : GEN COND, t1 = gsubgs(t, 1), BAD = primedivisors(hgm_get_BAD(hgm));
1755 203 : long k = hgm_get_WT(hgm) + 1, i, j, l;
1756 :
1757 203 : if (isintzero(t1)) T1 = cgetg(1, t_VECSMALL);
1758 : else
1759 : {
1760 196 : GEN fa = absZ_factor(numer_i(t1)), P = gel(fa,1), E = gel(fa,2);
1761 196 : if (!odd(k)) T1 = removebad(ZV_to_zv(P), BAD);
1762 : else
1763 : {
1764 182 : long j, lP = lg(P);
1765 182 : T1 = cgetg(lP, t_VECSMALL);
1766 350 : for (i = j = 1; i < lP; i++)
1767 : {
1768 168 : long p = itos(gel(P,i));
1769 168 : if (mpodd(gel(E,i)) && !zv_search(BAD, p)) T1[j++] = p;
1770 : }
1771 182 : setlg(T1,j);
1772 : }
1773 : }
1774 203 : COND = zv_prod_Z(T1);
1775 203 : T0 = removebad(gconcat(primedivisors(numer_i(t)),
1776 : primedivisors(denom_i(t))), BAD);
1777 217 : for (i = 1; i < lg(T0); i++)
1778 : {
1779 14 : long p = T0[i], e = tameexpo(hgm, Q_lval(t, p));
1780 14 : COND = mulii(powuu(p, e), COND);
1781 : }
1782 203 : l = lg(BAD); v = cgetg(l, t_VEC);
1783 490 : for (i = j = 1; i < l; i++) /* FIXME: precompute wild Euler factors */
1784 287 : if (hgmclass(hgm, BAD[i], t) == C_BAD)
1785 252 : gel(v,j++) = mkvec2(utoipos(BAD[i]), pol_1(0));
1786 203 : setlg(v, j); return gerepilecopy(av, mkvec4(VGA, utoi(k), v, COND));
1787 : }
1788 :
1789 : /***************************************************************/
1790 : /* PART IX: GUESS OF WILD PRIMES */
1791 : /***************************************************************/
1792 : static GEN
1793 1708 : vecmellin(GEN L, GEN K, GEN t, GEN td, GEN vroots, long bitprec)
1794 : {
1795 1708 : long i, N = lfunthetacost(L, t, 0, bitprec);
1796 1708 : GEN v = cgetg(N+1, t_VEC);
1797 427231 : for (i = 1; i <= N; i++)
1798 425523 : gel(v,i) = gammamellininvrt(K, gmul(td, gel(vroots,i)), bitprec);
1799 1708 : return v;
1800 : }
1801 :
1802 : /* List of Weil polynomials for prime p of degree d and weight w */
1803 : static GEN
1804 1610 : listweil_i(ulong d, ulong p, ulong w)
1805 : {
1806 : GEN V;
1807 1610 : if (d == 0) return mkvec(pol_1(0));
1808 1225 : if (odd(d))
1809 : {
1810 : GEN P;
1811 644 : if (odd(w)) return cgetg(1, t_VEC);
1812 420 : V = listweil_i(d - 1, p, w);
1813 420 : P = monomial(powuu(p, w / 2), 1, 0);
1814 420 : return shallowconcat(gmul(gsubsg(1, P), V), gmul(gaddsg(1, P), V));
1815 : }
1816 581 : if (d == 2)
1817 : {
1818 553 : long q = upowuu(p, w), s = usqrt(4 * q), i;
1819 553 : GEN Q = utoi(q);
1820 553 : V = cgetg(2*s + 3, t_VEC);
1821 4844 : for (i = 1; i <= 2*s + 1; i++) gel(V,i) = mkpoln(3, Q, stoi(1+s-i), gen_1);
1822 553 : gel(V, 2*s + 2) = mkpoln(3, negi(Q), gen_0, gen_1);
1823 553 : return V;
1824 : }
1825 28 : if (d == 4)
1826 : {
1827 28 : long q = upowuu(p, w), s = usqrt(16 * q), a, c, is2 = usqrt(4 * q);
1828 28 : double s2 = 2 * sqrt((double)q);
1829 28 : GEN Q2 = sqru(q), W, A, mA, tmp, Q;
1830 28 : V = cgetg(s + 3, t_VEC);
1831 252 : for (a = 0; a <= s; a++)
1832 : {
1833 224 : long b, m = ceil(a * s2) - 2 * q, M = ((a * a) >> 2) + 2 * q;
1834 224 : GEN AQ = stoi(a * q), mAQ = stoi(-a * q);
1835 224 : A = stoi(a); mA = stoi(-a);
1836 224 : W = cgetg(2*(M - m) + 3, t_VEC);
1837 1876 : for (b = m, c = 1; b <= M; b++)
1838 : {
1839 1652 : if (a) gel(W, c++) = mkpoln(5, Q2, mAQ, stoi(b), mA, gen_1);
1840 1652 : gel(W, c++) = mkpoln(5, Q2, AQ, stoi(b), A, gen_1);
1841 : }
1842 224 : setlg(W, c); gel(V, a + 1) = W;
1843 : }
1844 28 : W = cgetg(2 * is2 + 2, t_VEC);
1845 28 : tmp = mkpoln(3, stoi(-q), gen_0, gen_1);
1846 28 : Q = utoipos(q);
1847 147 : for (a = 0, c = 1; a <= is2; a++)
1848 : {
1849 119 : A = stoi(a); mA = stoi(-a);
1850 119 : if (a) gel(W, c++) = gmul(tmp, mkpoln(3, Q, mA, gen_1));
1851 119 : gel(W, c++) = gmul(tmp, mkpoln(3, Q, A, gen_1));
1852 : }
1853 28 : setlg(W, c); gel(V, s + 2) = W;
1854 28 : return shallowconcat1(V);
1855 : }
1856 0 : pari_err_IMPL("d > 5 in listweil");
1857 : return NULL; /* LCOV_EXCL_LINE */
1858 : }
1859 :
1860 : /* Same, weight <= w, by decreasing weight */
1861 : static GEN
1862 490 : listweilallw_i(ulong d, ulong p, ulong w)
1863 : {
1864 490 : long i = d? w: 0;
1865 490 : GEN V = cgetg(i+2, t_VEC);
1866 1680 : for (; i >= 0; i--) gel(V,i+1) = listweil_i(d, p, i);
1867 490 : return shallowconcat1(V);
1868 : }
1869 :
1870 : static long
1871 127344 : sumdeg(GEN W, GEN M)
1872 : {
1873 127344 : long i, l = lg(M), s = 0;
1874 375851 : for (i = 1; i < l; i++) s += degpol(gmael(W,i,M[i]+1));
1875 127344 : return s;
1876 : }
1877 :
1878 : static GEN
1879 35 : mul(GEN a, GEN b, ulong lim)
1880 : {
1881 35 : long na = lg(a)-1, nb = lg(b)-1, i, j, n;
1882 35 : GEN c = cgetg(na * nb + 1, t_VECSMALL);
1883 378 : for (i = n = 1; i <= na; i++)
1884 2597 : for (j = 1; j <= nb; j++)
1885 : {
1886 2254 : ulong m = umuluu_or_0(a[i], b[j]);
1887 2254 : if (m && m <= lim) c[n++] = m;
1888 : }
1889 35 : setlg(c, n); return c;
1890 : }
1891 : static GEN
1892 98 : listcond(GEN BAD, GEN achi, ulong min, ulong max)
1893 : {
1894 98 : long i, j, l = lg(BAD);
1895 98 : GEN P = cgetg(l, t_VEC), z;
1896 231 : for (i = 1; i < l; i++)
1897 : {
1898 133 : long p = BAD[i], v = achi[i];
1899 133 : gel(P,i) = upowers_u(p, ulogint(max, p) - v, upowuu(p, v));
1900 : }
1901 133 : z = gel(P,1); for (i = 2; i < l; i++) z = mul(z, gel(P,i), max);
1902 98 : if (min > 1)
1903 : {
1904 0 : l = lg(z);
1905 0 : for (i = j = 1; i < l; i++)
1906 0 : if ((ulong)z[i] >= min) z[j++] = z[i];
1907 0 : setlg(z, j);
1908 : }
1909 98 : vecsmall_sort(z); return z;
1910 : }
1911 :
1912 : /* Artificial lfundiv by zeta(s - (k - 1) / 2). */
1913 : static GEN
1914 35 : lfundivraw(GEN L)
1915 : {
1916 35 : long k = itos(ldata_get_k(L)), i;
1917 35 : GEN v = ldata_get_gammavec(L);
1918 35 : i = RgV_isin(v, stoi(-(k - 1) >> 2));
1919 35 : if (!i) pari_err_BUG("lfundivraw");
1920 35 : L = shallowcopy(L);
1921 35 : gel(L, 3) = vecsplice(v, i);
1922 35 : setlg(L, 7); return L;
1923 : }
1924 :
1925 : /* Compute forvec vectors from v, sorted by increasing total degree */
1926 : static GEN
1927 98 : forvecsort(GEN vF, GEN v)
1928 : {
1929 98 : GEN w, E = cyc2elts(v);
1930 98 : long i, l = lg(E);
1931 98 : w = cgetg(l, t_VECSMALL);
1932 127442 : for (i = 1; i < l; i++) w[i] = sumdeg(vF, gel(E, i));
1933 98 : return vecpermute(E, vecsmall_indexsort(w)); /* by increasing total degree */
1934 : }
1935 :
1936 : static GEN
1937 98 : lfunhgmwild(GEN L, GEN H, GEN t, GEN BAD, long pole, GEN hint, long bitprec)
1938 : {
1939 : GEN v, K, t0, t0r, t0ir, t0i, t0k, N0, vM, vD, val, PPOL, vF, achi;
1940 98 : long d, lM, iN, iM, i, k, k2, prec = nbits2prec(bitprec), lB = lg(BAD);
1941 : long BADprod, limdeg;
1942 98 : ulong minN = 1, maxN = 2048;
1943 :
1944 98 : v = cgetg(lB, t_VECSMALL); PPOL = cgetg(lB, t_VEC);
1945 231 : for (i = 1; i < lB; i++)
1946 : {
1947 133 : v[i] = itou( gmael(BAD,i,1) );
1948 133 : gel(PPOL,i) = shallowcopy(gel(BAD,i));
1949 : }
1950 98 : BAD = v;
1951 98 : BADprod = zv_prod(BAD);
1952 98 : achi = get_achi(H, t, BAD);
1953 98 : if (pole) L = lfundivraw(L);
1954 98 : limdeg = d = ldata_get_degree(L);
1955 98 : N0 = ldata_get_conductor(L);
1956 98 : if (hint) switch(typ(hint))
1957 : {
1958 : long l;
1959 0 : case t_INT:
1960 0 : limdeg = itos(hint);
1961 0 : if (limdeg < 0 || limdeg > d) pari_err_TYPE("lfunhgm [bad hint]", hint);
1962 0 : break;
1963 0 : case t_VEC:
1964 0 : l = lg(hint);
1965 0 : if (l > 1 && l < 4)
1966 : {
1967 0 : GEN t = gel(hint,1), r;
1968 0 : if (typ(t) != t_INT || signe(t) <= 0)
1969 0 : pari_err_TYPE("lfunhgm [bad hint]", hint);
1970 0 : t = dvmdii(t, N0, &r);
1971 0 : if (r != gen_0)
1972 0 : pari_err_TYPE("lfunhgm [bad hint]", hint);
1973 0 : minN = maxN = itou(t);
1974 0 : if (l == 3)
1975 : {
1976 0 : t = gel(hint,2);
1977 0 : if (typ(t) != t_INT || signe(t) < 0 || cmpis(t, d) > 0)
1978 0 : pari_err_TYPE("lfunhgm [bad hint]", hint);
1979 0 : limdeg = itos(t);
1980 : }
1981 : }
1982 : }
1983 98 : k = itos(ldata_get_k(L)); k2 = (k-1) >> 1;
1984 98 : K = gammamellininvinit(ldata_get_gammavec(L), 0, bitprec + 32);
1985 98 : t0 = sstoQ(11, 10); t0i = ginv(t0); t0k = gpowgs(t0, k);
1986 98 : t0r = gpow(t0, sstoQ(2,d), prec); t0ir = ginv(t0r);
1987 : /* vF[i] list of Weil polynomials F for p = BAD[i],
1988 : * F = prod_j Fj, |reciprocal root of Fj|^2 = p^(w + 1 - nj)
1989 : * 2v_p(lc(F)) = deg F * (w + 1) - sum_j deg Fj * nj; */
1990 98 : vF = cgetg(lB, t_VEC);
1991 98 : vD = cgetg(lB, t_VEC); /* vD[k][l] = sum_j deg Fj * nj for F = vF[k][l] */
1992 98 : v = cgetg(lB, t_VECSMALL);
1993 231 : for (i = 1; i < lB; i++)
1994 : {
1995 133 : GEN W = cgetg(limdeg+2, t_VEC), D;
1996 133 : long p = BAD[i], j, lW;
1997 623 : for (j = 0; j <= limdeg; j++) gel(W,j+1) = listweilallw_i(j, p, k-1);
1998 133 : gel(vF, i) = W = shallowconcat1(W);
1999 133 : lW = lg(W); v[i] = lW-1;
2000 133 : gel(vD,i) = D = cgetg(lW, t_VEC);
2001 10136 : for (j = 1; j < lW; j++)
2002 : {
2003 10003 : GEN F = gel(W,j);
2004 10003 : D[j] = degpol(F) * k - 2 * Z_lval(leading_coeff(F), p);
2005 : }
2006 : }
2007 98 : vM = forvecsort(vF, v); lM = lg(vM);
2008 98 : if (DEBUGLEVEL) { err_printf(" lM = %ld ", lM); err_flush(); }
2009 98 : L = shallowcopy(L);
2010 98 : val = cgetg(lB, t_VECSMALL);
2011 : for(;;)
2012 0 : {
2013 : GEN z, vroots, an0, vN;
2014 : long lN, lim;
2015 :
2016 98 : z = listcond(BAD, achi, minN, maxN);
2017 98 : if (maxN == minN) /* from hint */
2018 : {
2019 0 : minN = 1; /* in case it fails */
2020 0 : maxN--;
2021 : }
2022 : else
2023 : {
2024 98 : minN = maxN+1;
2025 98 : maxN <<= 2;
2026 : }
2027 98 : lN = lg(z);
2028 98 : if (lN == 1) continue;
2029 98 : vN = cgetg(lN, t_VEC);
2030 1785 : for (i = 1; i < lN; i++) gel(vN, i) = muliu(N0, z[i]);
2031 98 : gel(L, 5) = gel(vN, lN - 1);
2032 98 : if (DEBUGLEVEL) err_printf("\nmaxN = %ld\n", itos(gel(L,5)));
2033 98 : lim = lfunthetacost(L, t0i, 0, bitprec);
2034 98 : vroots = dirpowers(lim, sstoQ(2, d), prec + EXTRAPREC64);
2035 98 : an0 = hgmcoefs(H, t, lim);
2036 98 : if (pole)
2037 : {
2038 35 : GEN w = vecpowuu(lim, k2);
2039 44450 : for (i = 1; i <= lim; i++)
2040 44415 : if (cgcd(i, BADprod) > 1) gel(w, i) = gen_0;
2041 35 : an0 = dirdiv(an0, w);
2042 : }
2043 854 : for (iN = 1; iN < lN; iN++)
2044 : {
2045 854 : pari_sp av0 = avma, av2;
2046 854 : GEN vecK, vecKi, N = gel(vN,iN), isqN = gpow(N, sstoQ(-1,d), prec);
2047 854 : if (DEBUGLEVEL) { err_printf(" %ld ", itos(N)); err_flush(); }
2048 854 : gel(L,5) = N;
2049 854 : vecK = vecmellin(L, K, t0, gmul(t0r, isqN), vroots, bitprec);
2050 854 : vecKi= vecmellin(L, K, t0i,gmul(t0ir, isqN), vroots, bitprec);
2051 2310 : for (i = 1; i < lB; i++) val[i] = Z_lval(N, BAD[i]);
2052 854 : setlg(an0, lg(vecKi)); av2 = avma;
2053 2305926 : for (iM = 1; iM < lM; iM++, set_avma(av2))
2054 : {
2055 2305170 : GEN M = gel(vM, iM), an, eno, den;
2056 3214589 : for (i = 1; i < lB; i++)
2057 : {
2058 2986858 : GEN F = gmael(vF, i, M[i]+1);
2059 2986858 : long D, dF = degpol(F), a = val[i];
2060 2986858 : if (a < d - dF) break;
2061 2797417 : if (a)
2062 : {
2063 2606261 : if (d == dF) break;
2064 1163771 : D = mael(vD, i, M[i]+1); /* sum_j deg Fj nj */
2065 1163771 : if (D == d && d - dF != a) break;
2066 : /* a = a(chi) + d - deg F - 1 */
2067 1100575 : if (D == d-1 && a != d - dF + achi[i] - 1 && !gequal1(t)) break;
2068 : }
2069 909419 : gmael(PPOL, i, 2) = F;
2070 : }
2071 2305170 : if (i < lB) continue;
2072 227731 : an = hgmmodif(an0, PPOL);
2073 227731 : den = gmul(t0k, RgV_dotproduct(vecK, an));
2074 227731 : eno = gdiv(RgV_dotproduct(vecKi, an), den);
2075 227731 : if (gexpo(gsubgs(gabs(eno, prec), 1)) < -bitprec / 2)
2076 : {
2077 98 : if (pole)
2078 77 : for (i = 1; i < lB; i++)
2079 : {
2080 42 : GEN t = deg1pol_shallow(negi(powuu(BAD[i], k2)), gen_1, 0);
2081 42 : gmael(PPOL, i, 2) = gmul(gmael(PPOL, i, 2), t);
2082 : }
2083 231 : for (i = 1; i < lB; i++) gmael(PPOL, i, 2) = ginv(gmael(PPOL, i, 2));
2084 98 : eno = grndtoi(eno, NULL);
2085 98 : if (typ(eno) != t_INT) pari_err_BUG("lfunhgmwild");
2086 98 : return mkvec3(N, mkvec2(t, PPOL), eno);
2087 : }
2088 : }
2089 756 : set_avma(av0);
2090 : }
2091 : }
2092 : }
2093 :
2094 : static GEN
2095 91 : BAD2small(GEN BAD)
2096 : {
2097 91 : long i, l = lg(BAD);
2098 91 : GEN v = cgetg(l, t_VECSMALL);
2099 210 : for (i = 1; i < l; i++) v[i] = itou( gmael(BAD,i,1) );
2100 91 : return v;
2101 : }
2102 :
2103 : /* moments of motive */
2104 : static GEN
2105 91 : hgmmoments(GEN H, GEN t, GEN M, long nb)
2106 : {
2107 91 : pari_sp av = avma, av2;
2108 91 : GEN P, v, S, L = hgmlfuninfty(H, t), BAD = BAD2small(gel(L, 3));
2109 91 : GEN k2 = gmul2n(gsubgs(gel(L,2), 1), -1);
2110 91 : long ct = 0, i, j, lP, lm, tm = typ(M);
2111 :
2112 91 : if (!nb) nb = 1000;
2113 91 : P = primes_zv(nb); lP = lg(P);
2114 91 : v = hgmcoefs(H, t, P[lP - 1]);
2115 91 : if (tm != t_VEC) M = gtovec(M);
2116 91 : av2 = avma; lm = lg(M);
2117 91 : S = const_vec(lm - 1, real_0(DEFAULTPREC));
2118 3731 : for (i = 1; i < lP; i++)
2119 : {
2120 3640 : long p = P[i];
2121 3640 : if (!Q_lval(t, p) && !zv_search(BAD, p))
2122 : {
2123 3500 : GEN T = gdiv(gel(v, p), gpow(utoipos(p), k2, DEFAULTPREC));
2124 3500 : ct++;
2125 7000 : for (j = 1; j < lm; j++)
2126 3500 : gel(S,j) = gadd(gel(S,j), gpow(T, gel(M,j), DEFAULTPREC));
2127 : }
2128 3640 : if ((i & 0xf) == 0) S = gerepilecopy(av2, S);
2129 : }
2130 91 : if (tm != t_VEC && tm != t_COL && tm != t_VECSMALL) S = gel(S, 1);
2131 91 : return gerepileupto(av, gdivgu(S, ct));
2132 : }
2133 :
2134 : /* Heuristic guess: is there a pole ? */
2135 : static long
2136 91 : lfunhgmispole(GEN H, GEN t, long nb)
2137 : {
2138 91 : pari_sp av = avma;
2139 : GEN M;
2140 91 : if (!nb) nb = 40;
2141 91 : M = hgmmoments(H, t, gen_1, nb);
2142 91 : return gc_bool(av, gexpo(M) > -2);
2143 : }
2144 :
2145 : static GEN
2146 112 : tag(GEN x, long t) { return mkvec2(mkvecsmall(t), x); }
2147 :
2148 : static GEN
2149 112 : lfunhgm_i(GEN hgm, GEN t, GEN hint, long bitprec)
2150 : {
2151 112 : GEN L, vr, v = hgmlfuninfty(hgm, t), vga = zv_to_ZV(gel(v,1)), k = gel(v,2);
2152 112 : GEN BAD = gel(v,3), COND = gel(v, 4);
2153 112 : long pole = mpodd(k) && lfunhgmispole(hgm, t, 0), lB = lg(BAD);
2154 :
2155 112 : L = mkvecn(pole? 7: 6, tag(mkvec2(hgm,t), t_LFUN_HGM), gen_0, vga, k, COND,
2156 : gen_0, mkvec(mkvec2(shifti(addiu(k,1),-1), gen_0)));
2157 112 : if (pole && ldata_get_degree(L) == 1)
2158 : { /* motive = zeta */
2159 : long i;
2160 0 : gel(L,5) = gel(L,6) = gel(L,7) = gen_1;
2161 0 : if (lB != 1)
2162 : {
2163 0 : GEN R = mkrfrac(gen_1, deg1pol_shallow(gen_m1,gen_1,0));
2164 0 : gmael3(L, 1, 2, 2) = mkvec2(t, BAD);
2165 0 : for (i = 1; i < lB; i++) gmael(BAD, i, 2) = R;
2166 : }
2167 0 : return L;
2168 : }
2169 112 : if (lB == 1)
2170 : {
2171 14 : vr = lfunrootres(L, bitprec);
2172 14 : gel(L, 6) = gel(vr, 3);
2173 14 : if (pole) gel(L, 7) = gel(vr, 1);
2174 14 : return L;
2175 : }
2176 98 : v = lfunhgmwild(L, hgm, t, BAD, pole, hint, bitprec);
2177 98 : gel(L, 5) = gel(v, 1); /* N */
2178 98 : gmael3(L, 1, 2, 2) = gel(v, 2); /* [t, PPOL] */
2179 98 : gel(L, 6) = gel(v, 3); /* w */
2180 98 : if (pole)
2181 : {
2182 35 : vr = lfunrootres(L, bitprec);
2183 35 : gel(L, 7) = gel(vr, 1);
2184 : }
2185 98 : return L;
2186 : }
2187 : GEN
2188 119 : lfunhgm(GEN hgm, GEN t, GEN hint, long bit)
2189 : {
2190 119 : pari_sp av = avma;
2191 119 : if (!checkhgm(hgm)) pari_err_TYPE("lfunhgm", hgm);
2192 112 : return gerepilecopy(av, lfunhgm_i(hgm, t, hint, bit));
2193 : }
2194 :
2195 : GEN
2196 3997 : dirhgm_worker(GEN P, ulong X, GEN hgm, GEN t)
2197 : {
2198 3997 : pari_sp av = avma;
2199 3997 : long i, l = lg(P);
2200 3997 : GEN W = cgetg(l, t_VEC);
2201 40822 : for(i = 1; i < l; i++)
2202 : {
2203 36825 : ulong p = uel(P,i);
2204 36825 : long e, d = ulogint(X, p) + 1; /* minimal d such that p^d > X */
2205 36825 : gel(W,i) = RgXn_inv(hgmeulerfactorlimit(hgm, t, p, d-1, 0, &e), d);
2206 : }
2207 3997 : return gerepilecopy(av, mkvec2(P,W));
2208 : }
2209 :
2210 : GEN
2211 399 : hgmcoefs(GEN hgm, GEN t, long n)
2212 : {
2213 399 : GEN worker, bad = NULL;
2214 399 : if (!checkhgm(hgm)) pari_err_TYPE("hgmcoefs", hgm);
2215 399 : if (typ(t) == t_VEC && lg(t) == 3) { bad = gel(t,2); t = gel(t,1); }
2216 399 : if (!is_rational_t(typ(t))) pari_err_TYPE("hgmcoefs",t);
2217 392 : worker = snm_closure(is_entry("_dirhgm_worker"), mkvec2(hgm, t));
2218 392 : return pardireuler(worker, gen_2, stoi(n), NULL, bad);
2219 : }
|
