Line data Source code
1 : /* Copyright (C) 2015 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; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : /********************************************************************/
16 : /** **/
17 : /** L-functions: Applications **/
18 : /** **/
19 : /********************************************************************/
20 :
21 : #include "pari.h"
22 : #include "paripriv.h"
23 :
24 : #define DEBUGLEVEL DEBUGLEVEL_lfun
25 :
26 : static GEN
27 26047 : tag(GEN x, long t) { return mkvec2(mkvecsmall(t), x); }
28 :
29 : /* v a t_VEC of length > 1 */
30 : static int
31 104971 : is_tagged(GEN v)
32 : {
33 104971 : GEN T = gel(v,1);
34 104971 : return (typ(T)==t_VEC && lg(T)==3 && typ(gel(T,1))==t_VECSMALL);
35 : }
36 : /* rough check */
37 : static long
38 126440 : is_ldata(GEN L)
39 : {
40 126440 : long l = lg(L);
41 126440 : return typ(L) == t_VEC && (l == 7 || l == 8);
42 : }
43 : /* thorough check */
44 : static void
45 104859 : checkldata(GEN ldata)
46 : {
47 : GEN vga, w, N;
48 : #if 0 /* assumed already checked and true */
49 : if (!is_ldata(ldata) || !is_tagged(ldata)) pari_err_TYPE("checkldata", ldata);
50 : #endif
51 104859 : vga = ldata_get_gammavec(ldata);
52 104859 : if (typ(vga) != t_VEC) pari_err_TYPE("checkldata [gammavec]",vga);
53 104859 : w = gel(ldata, 4); /* FIXME */
54 104859 : switch(typ(w))
55 : {
56 102829 : case t_INT: case t_FRAC: break;
57 2030 : case t_VEC: if (lg(w) == 3 && is_rational_t(typ(gel(w,1)))) break;
58 0 : default: pari_err_TYPE("checkldata [weight]",w);
59 : }
60 104859 : N = ldata_get_conductor(ldata);
61 104859 : if (typ(N) != t_INT) pari_err_TYPE("checkldata [conductor]",N);
62 104859 : }
63 :
64 : /* tag as t_LFUN_GENERIC */
65 : static void
66 609 : lfuncreate_tag(GEN L)
67 : {
68 609 : if (is_tagged(L)) return;
69 448 : gel(L,1) = tag(gel(L,1), t_LFUN_GENERIC);
70 448 : if (typ(gel(L,2)) != t_INT) gel(L,2) = tag(gel(L,2), t_LFUN_GENERIC);
71 : }
72 :
73 : /* shallow */
74 : static GEN
75 168 : closure2ldata(GEN C, long prec)
76 : {
77 168 : GEN L = closure_callgen0prec(C, prec);
78 168 : if (is_ldata(L)) { checkldata(L); lfuncreate_tag(L); }
79 56 : else L = lfunmisc_to_ldata_shallow(L);
80 168 : return L;
81 : }
82 :
83 : /* data may be either an object (polynomial, elliptic curve, etc...)
84 : * or a description vector [an,sd,Vga,k,conductor,rootno,{poles}]. */
85 : GEN
86 4025 : lfuncreate(GEN data)
87 : {
88 4025 : if (is_ldata(data))
89 : {
90 497 : GEN L = gcopy(data);
91 497 : lfuncreate_tag(L); checkldata(L); return L;
92 : }
93 3528 : if (typ(data) == t_CLOSURE && closure_arity(data)==0)
94 : {
95 14 : pari_sp av = avma;
96 14 : GEN L = closure2ldata(data, DEFAULTPREC);
97 14 : gel(L,1) = tag(data, t_LFUN_CLOSURE0); return gerepilecopy(av, L);
98 : }
99 3514 : return lfunmisc_to_ldata(data);
100 : }
101 :
102 : GEN
103 133 : lfunparams(GEN L, long prec)
104 : {
105 133 : pari_sp av = avma;
106 : GEN k, N, v;
107 : long p;
108 :
109 133 : if (!is_ldata(L) || !is_tagged(L)) L = lfunmisc_to_ldata_shallow(L);
110 133 : N = ldata_get_conductor(L);
111 133 : k = ldata_get_k(L);
112 133 : v = ldata_get_gammavec(L);
113 133 : p = gprecision(v);
114 133 : if (p > prec) v = gprec_wtrunc(v, prec);
115 133 : else if (p < prec)
116 : {
117 133 : GEN van = ldata_get_an(L), an = gel(van,2);
118 133 : long t = mael(van,1,1);
119 133 : if (t == t_LFUN_CLOSURE0) L = closure2ldata(an, prec);
120 : }
121 133 : return gerepilecopy(av, mkvec3(N, k, v));
122 : }
123 :
124 : /********************************************************************/
125 : /** Simple constructors **/
126 : /********************************************************************/
127 : static GEN ldata_eulerf(GEN van, GEN p, long prec);
128 :
129 : static GEN
130 126 : vecan_conj(GEN an, long n, long prec)
131 : {
132 126 : GEN p1 = ldata_vecan(gel(an,1), n, prec);
133 126 : return typ(p1) == t_VEC? conj_i(p1): p1;
134 : }
135 :
136 : static GEN
137 0 : eulerf_conj(GEN an, GEN p, long prec)
138 : {
139 0 : GEN p1 = ldata_eulerf(gel(an,1), p, prec);
140 0 : return conj_i(p1);
141 : }
142 :
143 : static GEN
144 308 : vecan_mul(GEN an, long n, long prec)
145 : {
146 308 : GEN p1 = ldata_vecan(gel(an,1), n, prec);
147 308 : GEN p2 = ldata_vecan(gel(an,2), n, prec);
148 308 : if (typ(p1) == t_VECSMALL) p1 = vecsmall_to_vec(p1);
149 308 : if (typ(p2) == t_VECSMALL) p2 = vecsmall_to_vec(p2);
150 308 : return dirmul(p1, p2);
151 : }
152 :
153 : static GEN
154 28 : eulerf_mul(GEN an, GEN p, long prec)
155 : {
156 28 : GEN p1 = ldata_eulerf(gel(an,1), p, prec);
157 28 : GEN p2 = ldata_eulerf(gel(an,2), p, prec);
158 28 : return gmul(p1, p2);
159 : }
160 :
161 : static GEN
162 77 : lfunconvol(GEN a1, GEN a2)
163 77 : { return tag(mkvec2(a1, a2), t_LFUN_MUL); }
164 :
165 : static GEN
166 630 : vecan_div(GEN an, long n, long prec)
167 : {
168 630 : GEN p1 = ldata_vecan(gel(an,1), n, prec);
169 630 : GEN p2 = ldata_vecan(gel(an,2), n, prec);
170 630 : if (typ(p1) == t_VECSMALL) p1 = vecsmall_to_vec(p1);
171 630 : if (typ(p2) == t_VECSMALL) p2 = vecsmall_to_vec(p2);
172 630 : return dirdiv(p1, p2);
173 : }
174 :
175 : static GEN
176 14 : eulerf_div(GEN an, GEN p, long prec)
177 : {
178 14 : GEN p1 = ldata_eulerf(gel(an,1), p, prec);
179 14 : GEN p2 = ldata_eulerf(gel(an,2), p, prec);
180 14 : return gdiv(p1, p2);
181 : }
182 :
183 : static GEN
184 56 : lfunconvolinv(GEN a1, GEN a2)
185 56 : { return tag(mkvec2(a1,a2), t_LFUN_DIV); }
186 :
187 : static GEN
188 84 : lfunconj(GEN a1)
189 84 : { return tag(mkvec(a1), t_LFUN_CONJ); }
190 :
191 : static GEN
192 133 : lfuncombdual(GEN (*fun)(GEN, GEN), GEN ldata1, GEN ldata2)
193 : {
194 133 : GEN a1 = ldata_get_an(ldata1), a2 = ldata_get_an(ldata2);
195 133 : GEN b1 = ldata_get_dual(ldata1), b2 = ldata_get_dual(ldata2);
196 133 : if (typ(b1)==t_INT && typ(b2)==t_INT)
197 133 : return utoi(signe(b1) || signe(b2));
198 : else
199 : {
200 0 : if (typ(b1)==t_INT) b1 = signe(b1) ? lfunconj(a1): a1;
201 0 : if (typ(b2)==t_INT) b2 = signe(b2) ? lfunconj(a2): a2;
202 0 : return fun(b1, b2);
203 : }
204 : }
205 :
206 : static GEN
207 2877 : vecan_twist(GEN an, long n, long prec)
208 : {
209 2877 : GEN p1 = ldata_vecan(gel(an,1), n, prec);
210 2877 : GEN p2 = ldata_vecan(gel(an,2), n, prec);
211 : long i;
212 : GEN V;
213 2877 : if (typ(p1) == t_VECSMALL) p1 = vecsmall_to_vec(p1);
214 2877 : if (typ(p2) == t_VECSMALL) p2 = vecsmall_to_vec(p2);
215 2877 : V = cgetg(n+1, t_VEC);
216 1424479 : for(i = 1; i <= n ; i++)
217 1421602 : gel(V, i) = gmul(gel(p1, i), gel(p2, i));
218 2877 : return V;
219 : }
220 :
221 : static GEN
222 14 : eulerf_twist(GEN an, GEN p, long prec)
223 : {
224 14 : GEN p1 = ldata_eulerf(gel(an,1), p, prec);
225 14 : GEN p2 = ginv(ldata_eulerf(gel(an,2), p, prec));
226 14 : if (typ(p2)!=t_POL || degpol(p2)==0)
227 0 : return poleval(p1,pol_0(0));
228 14 : if (degpol(p2)!=1) pari_err_IMPL("lfuneuler");
229 14 : return poleval(p1,monomial(gneg(gel(p2,3)),1,0));
230 : }
231 :
232 : static GEN
233 637 : vecan_shift(GEN an, long n, long prec)
234 : {
235 637 : GEN p1 = ldata_vecan(gel(an,1), n, prec);
236 637 : GEN s = gel(an,2);
237 : long i;
238 : GEN V;
239 637 : if (typ(p1) == t_VECSMALL) p1 = vecsmall_to_vec(p1);
240 637 : V = cgetg(n+1, t_VEC);
241 637 : if (typ(s)==t_INT)
242 : {
243 448 : if (equali1(s))
244 31850 : for(i = 1; i <= n ; i++)
245 : {
246 31500 : GEN gi = gel(p1, i);
247 31500 : gel(V, i) = gequal0(gi)? gi: gmulgu(gi, i);
248 : }
249 : else
250 1848 : for(i = 1; i <= n ; i++)
251 : {
252 1750 : GEN gi = gel(p1, i);
253 1750 : gel(V, i) = gequal0(gi)? gi: gmul(gi, powgi(utoi(i), s));
254 : }
255 : }
256 : else
257 : {
258 189 : GEN D = dirpowers(n, s, prec);
259 7049 : for(i = 1; i <= n ; i++)
260 6860 : gel(V, i) = gmul(gel(p1,i), gel(D,i));
261 : }
262 637 : return V;
263 : }
264 :
265 : static GEN
266 42 : eulerf_shift(GEN an, GEN p, long prec)
267 : {
268 42 : GEN p1 = ldata_eulerf(gel(an,1), p, prec);
269 42 : GEN s = gel(an,2);
270 42 : return gsubst(p1, 0, monomial(gpow(p, s, prec), 1, 0));
271 : }
272 :
273 : static GEN
274 84 : eulerf_hgm(GEN an, GEN p)
275 : {
276 84 : GEN H = gel(an,1), t = gel(an,2);
277 84 : if (typ(t)==t_VEC && lg(t)==3)
278 : {
279 28 : GEN L = gel(t,2);
280 28 : long i, l = lg(L);
281 28 : t = gel(t,1);
282 49 : for (i = 1; i < l; i++) /* wild primes */
283 35 : if (equalii(p, gmael(L, i, 1))) break;
284 28 : if (i<l) return gmael(L,i,2);
285 : }
286 70 : return ginv(hgmeulerfactor(H, t, itos(p), NULL));
287 : }
288 :
289 : static GEN
290 315 : deg1ser_shallow(GEN a1, GEN a0, long e)
291 315 : { return RgX_to_ser(deg1pol_shallow(a1, a0, 0), e+2); }
292 : /* lfunrtopoles without sort */
293 : static GEN
294 168 : rtopoles(GEN r)
295 : {
296 168 : long j, l = lg(r);
297 168 : GEN v = cgetg(l, t_VEC);
298 350 : for (j = 1; j < l; j++)
299 : {
300 182 : GEN rj = gel(r,j), a = gel(rj,1);
301 182 : gel(v,j) = a;
302 : }
303 168 : return v;
304 : }
305 : /* re = polar part; overestimate when re = gen_0 (unknown) */
306 : static long
307 266 : orderpole(GEN re) { return typ(re) == t_SER? -valser(re): 1; }
308 : static GEN
309 77 : lfunmulpoles(GEN ldata1, GEN ldata2, long bitprec)
310 : {
311 77 : GEN r, k = ldata_get_k(ldata1), b1 = NULL, b2 = NULL;
312 77 : GEN r1 = ldata_get_residue(ldata1);
313 77 : GEN r2 = ldata_get_residue(ldata2);
314 77 : long i, j, l, L = 0;
315 :
316 77 : if (!r1 && !r2) return NULL;
317 70 : if (r1 && !is_vec_t(typ(r1))) r1 = mkvec(mkvec2(k, r1));
318 70 : if (r2 && !is_vec_t(typ(r2))) r2 = mkvec(mkvec2(k, r2));
319 70 : if (r1) { b1 = rtopoles(r1); L += lg(b1); }
320 70 : if (r2) { b2 = rtopoles(r2); L += lg(b2); }
321 70 : r = cgetg(L, t_VEC); j = 1;
322 70 : if (b1)
323 : {
324 63 : l = lg(b1);
325 133 : for (i = 1; i < l; i++)
326 : {
327 70 : GEN z, z1, z2, be = gmael(r1,i,1);
328 70 : long n, v = orderpole(gmael(r1,i,2));
329 70 : if (b2 && (n = RgV_isin(b2, be))) v += orderpole(gmael(r2,n,2));
330 70 : z = deg1ser_shallow(gen_1, be, 2 + v);
331 70 : z1 = lfun(ldata1,z,bitprec);
332 70 : z2 = lfun(ldata2,z,bitprec);
333 70 : gel(r,j++) = mkvec2(be, gmul(z1, z2));
334 : }
335 : }
336 70 : if (b2)
337 : {
338 63 : long l = lg(b2);
339 133 : for (i = 1; i < l; i++)
340 : {
341 70 : GEN z, z1, z2, be = gmael(r2,i,1);
342 70 : long n, v = orderpole(gmael(r2,i,2));
343 70 : if (b1 && (n = RgV_isin(b1, be))) continue; /* done already */
344 28 : z = deg1ser_shallow(gen_1, be, 2 + v);
345 28 : z1 = lfun(ldata1,z,bitprec);
346 28 : z2 = lfun(ldata2,z,bitprec);
347 28 : gel(r,j++) = mkvec2(be, gmul(z1, z2));
348 : }
349 : }
350 70 : setlg(r, j); return r;
351 : }
352 :
353 : static GEN
354 77 : lfunmul_k(GEN ldata1, GEN ldata2, GEN k, long bitprec)
355 : {
356 : GEN r, N, Vga, eno, a1a2, b1b2;
357 77 : r = lfunmulpoles(ldata1, ldata2, bitprec);
358 77 : N = gmul(ldata_get_conductor(ldata1), ldata_get_conductor(ldata2));
359 77 : Vga = shallowconcat(ldata_get_gammavec(ldata1), ldata_get_gammavec(ldata2));
360 77 : Vga = sort(Vga);
361 77 : eno = gmul(ldata_get_rootno(ldata1), ldata_get_rootno(ldata2));
362 77 : a1a2 = lfunconvol(ldata_get_an(ldata1), ldata_get_an(ldata2));
363 77 : b1b2 = lfuncombdual(lfunconvol, ldata1, ldata2);
364 77 : return mkvecn(r? 7: 6, a1a2, b1b2, Vga, k, N, eno, r);
365 : }
366 :
367 : GEN
368 63 : lfunmul(GEN ldata1, GEN ldata2, long bitprec)
369 : {
370 63 : pari_sp ltop = avma;
371 : GEN k;
372 63 : long prec = nbits2prec(bitprec);
373 63 : ldata1 = ldata_newprec(lfunmisc_to_ldata_shallow(ldata1), prec);
374 63 : ldata2 = ldata_newprec(lfunmisc_to_ldata_shallow(ldata2), prec);
375 63 : k = ldata_get_k(ldata1);
376 63 : if (!gequal(ldata_get_k(ldata2),k))
377 0 : pari_err_OP("lfunmul [weight]",ldata1, ldata2);
378 63 : return gerepilecopy(ltop, lfunmul_k(ldata1, ldata2, k, bitprec));
379 : }
380 :
381 : static GEN
382 56 : lfundivpoles(GEN ldata1, GEN ldata2, long bitprec)
383 : {
384 : long i, j, l;
385 56 : GEN be2, k = ldata_get_k(ldata1);
386 56 : GEN r1 = ldata_get_residue(ldata1);
387 56 : GEN r2 = ldata_get_residue(ldata2), r;
388 :
389 56 : if (r1 && !is_vec_t(typ(r1))) r1 = mkvec(mkvec2(k, r1));
390 56 : if (r2 && !is_vec_t(typ(r2))) r2 = mkvec(mkvec2(k, r2));
391 56 : if (!r1) return NULL;
392 56 : l = lg(r1); r = cgetg(l, t_VEC);
393 56 : be2 = r2? rtopoles(r2): NULL;
394 112 : for (i = j = 1; j < l; j++)
395 : {
396 56 : GEN z, v = gel(r1,j), be = gel(v,1), s1 = gel(v,2);
397 : long n;
398 56 : if (be2 && (n = RgV_isin(be2, be)))
399 : {
400 42 : GEN s2 = gmael(r2,n,2); /* s1,s2: polar parts */
401 42 : if (orderpole(s1) == orderpole(s2)) continue;
402 : }
403 14 : z = gdiv(lfun(ldata1,be,bitprec), lfun(ldata2,be,bitprec));
404 14 : if (valser(z) < 0) gel(r,i++) = mkvec2(be, z);
405 : }
406 56 : if (i == 1) return NULL;
407 14 : setlg(r, i); return r;
408 : }
409 :
410 : static GEN
411 56 : lfunvgasub(GEN v01, GEN v2)
412 : {
413 56 : GEN v1 = shallowcopy(v01), v;
414 56 : long l1 = lg(v1), l2 = lg(v2), j1, j2, j;
415 119 : for (j2 = 1; j2 < l2; j2++)
416 : {
417 91 : for (j1 = 1; j1 < l1; j1++)
418 91 : if (gel(v1,j1) && gequal(gel(v1,j1), gel(v2,j2)))
419 : {
420 63 : gel(v1,j1) = NULL; break;
421 : }
422 63 : if (j1 == l1) pari_err_OP("lfunvgasub", v1, v2);
423 : }
424 56 : v = cgetg(l1-l2+1, t_VEC);
425 259 : for (j1 = j = 1; j1 < l1; j1++)
426 203 : if (gel(v1, j1)) gel(v,j++) = gel(v1,j1);
427 56 : return v;
428 : }
429 :
430 : GEN
431 56 : lfundiv(GEN ldata1, GEN ldata2, long bitprec)
432 : {
433 56 : pari_sp ltop = avma;
434 : GEN k, r, N, v, eno, a1a2, b1b2, eno2;
435 56 : long prec = nbits2prec(bitprec);
436 56 : ldata1 = ldata_newprec(lfunmisc_to_ldata_shallow(ldata1), prec);
437 56 : ldata2 = ldata_newprec(lfunmisc_to_ldata_shallow(ldata2), prec);
438 56 : k = ldata_get_k(ldata1);
439 56 : if (!gequal(ldata_get_k(ldata2),k))
440 0 : pari_err_OP("lfundiv [weight]",ldata1, ldata2);
441 56 : N = gdiv(ldata_get_conductor(ldata1), ldata_get_conductor(ldata2));
442 56 : if (typ(N) != t_INT) pari_err_OP("lfundiv [conductor]",ldata1, ldata2);
443 56 : r = lfundivpoles(ldata1, ldata2, bitprec);
444 56 : a1a2 = lfunconvolinv(ldata_get_an(ldata1), ldata_get_an(ldata2));
445 56 : b1b2 = lfuncombdual(lfunconvolinv, ldata1, ldata2);
446 56 : eno2 = ldata_get_rootno(ldata2);
447 56 : eno = isintzero(eno2)? gen_0: gdiv(ldata_get_rootno(ldata1), eno2);
448 56 : v = lfunvgasub(ldata_get_gammavec(ldata1), ldata_get_gammavec(ldata2));
449 56 : return gerepilecopy(ltop, mkvecn(r? 7: 6, a1a2, b1b2, v, k, N, eno, r));
450 : }
451 :
452 : static GEN
453 2387 : gamma_imagchi(GEN gam, GEN w)
454 : {
455 2387 : long i, j, k=1, l;
456 2387 : GEN g = cgetg_copy(gam, &l);
457 2387 : gam = shallowcopy(gam);
458 7161 : for (i = l-1; i>=1; i--)
459 : {
460 4774 : GEN al = gel(gam, i);
461 4774 : if (al)
462 : {
463 2394 : GEN N = gadd(w,gmul2n(real_i(al),1));
464 2394 : if (gcmpgs(N,2) > 0)
465 : {
466 2380 : GEN bl = gsubgs(al, 1);
467 2380 : for (j=1; j < i; j++)
468 2380 : if (gel(gam,j) && gequal(gel(gam,j), bl))
469 2380 : { gel(gam,j) = NULL; break; }
470 2380 : if (j==i) return NULL;
471 2380 : gel(g, k++) = al;
472 2380 : gel(g, k++) = bl;
473 14 : } else if (gequal0(N))
474 14 : gel(g, k++) = gaddgs(al, 1);
475 0 : else if (gequal1(N))
476 0 : gel(g, k++) = gsubgs(al, 1);
477 0 : else return NULL;
478 : }
479 : }
480 2387 : return sort(g);
481 : }
482 :
483 : GEN
484 5173 : lfuntwist(GEN ldata1, GEN chi, long bitprec)
485 : {
486 5173 : pari_sp ltop = avma;
487 : GEN k, L, N, N1, N2, a, a1, a2, b, b1, b2, gam, gam1, gam2;
488 : GEN ldata2;
489 : long d1, t;
490 5173 : long prec = nbits2prec(bitprec);
491 5173 : ldata1 = ldata_newprec(lfunmisc_to_ldata_shallow(ldata1), prec);
492 5173 : ldata2 = lfunmisc_to_ldata_shallow(chi);
493 5173 : t = ldata_get_type(ldata2);
494 5173 : a1 = ldata_get_an(ldata1);
495 5173 : a2 = ldata_get_an(ldata2);
496 5173 : if (t == t_LFUN_ZETA)
497 2415 : return gerepilecopy(ltop, ldata1);
498 2758 : if (t != t_LFUN_CHIZ && t != t_LFUN_KRONECKER &&
499 7 : ( t != t_LFUN_CHIGEN || nf_get_degree(bnr_get_nf(gmael(a2,2,1))) != 1))
500 0 : pari_err_TYPE("lfuntwist", chi);
501 2758 : N1 = ldata_get_conductor(ldata1);
502 2758 : N2 = ldata_get_conductor(ldata2);
503 2758 : if (!gequal1(gcdii(N1, N2)))
504 0 : pari_err_IMPL("lfuntwist (conductors not coprime)");
505 2758 : k = ldata_get_k(ldata1);
506 2758 : d1 = ldata_get_degree(ldata1);
507 2758 : N = gmul(N1, gpowgs(N2, d1));
508 2758 : gam1 = ldata_get_gammavec(ldata1);
509 2758 : gam2 = ldata_get_gammavec(ldata2);
510 2758 : if (gequal0(gel(gam2, 1)))
511 371 : gam = gam1;
512 : else
513 2387 : gam = gamma_imagchi(ldata_get_gammavec(ldata1), gaddgs(k,-1));
514 2758 : if (!gam) pari_err_IMPL("lfuntwist (gammafactors)");
515 2758 : b1 = ldata_get_dual(ldata1);
516 2758 : b2 = ldata_get_dual(ldata2);
517 2758 : a = tag(mkvec2(a1, a2), t_LFUN_TWIST);
518 2758 : if (typ(b1)==t_INT)
519 2758 : b = signe(b1) && signe(b2) ? gen_0: gen_1;
520 : else
521 0 : b = tag(mkvec2(b1,lfunconj(a2)), t_LFUN_TWIST);
522 2758 : L = mkvecn(6, a, b, gam, k, N, gen_0);
523 2758 : return gerepilecopy(ltop, L);
524 : }
525 :
526 : static GEN
527 280 : lfundualpoles(GEN ldata, GEN reno)
528 : {
529 : long l, j;
530 280 : GEN k = ldata_get_k(ldata);
531 280 : GEN r = gel(reno,2), eno = gel(reno,3), R;
532 280 : R = cgetg_copy(r, &l);
533 854 : for (j = 1; j < l; j++)
534 : {
535 574 : GEN b = gmael(r,j,1), e = gmael(r,j,2);
536 574 : long v = varn(e);
537 574 : GEN E = gsubst(gdiv(e, eno), v, gneg(pol_x(v)));
538 574 : gel(R,l-j) = mkvec2(gsub(k,b), E);
539 : }
540 280 : return R;
541 : }
542 :
543 : static GEN
544 567 : ginvvec(GEN x)
545 : {
546 567 : if (is_vec_t(typ(x)))
547 42 : pari_APPLY_same(ginv(gel(x,i)))
548 : else
549 553 : return ginv(x);
550 : }
551 :
552 : GEN
553 630 : lfundual(GEN L, long bitprec)
554 : {
555 630 : pari_sp av = avma;
556 630 : long prec = nbits2prec(bitprec);
557 630 : GEN ldata = ldata_newprec(lfunmisc_to_ldata_shallow(L), prec);
558 630 : GEN a = ldata_get_an(ldata), b = ldata_get_dual(ldata);
559 630 : GEN e = ldata_get_rootno(ldata);
560 630 : GEN ldual, ad, bd, ed, Rd = NULL;
561 630 : if (typ(b) == t_INT)
562 : {
563 602 : ad = equali1(b) ? lfunconj(a): a;
564 602 : bd = b;
565 : }
566 28 : else { ad = b; bd = a; }
567 630 : if (lg(ldata)==8)
568 : {
569 280 : GEN reno = lfunrootres(ldata, bitprec);
570 280 : e = gel(reno,3);
571 280 : Rd = lfundualpoles(ldata, reno);
572 : }
573 630 : ed = isintzero(e) ? e: ginvvec(e);
574 630 : ldual = mkvecn(Rd ? 7:6, ad, bd, gel(ldata,3), gel(ldata,4), gel(ldata,5), ed, Rd);
575 630 : return gerepilecopy(av, ldual);
576 : }
577 :
578 : static GEN
579 98 : RgV_Rg_translate(GEN x, GEN s)
580 259 : { pari_APPLY_same(gadd(gel(x,i),s)) }
581 :
582 : static GEN
583 28 : pole_translate(GEN x, GEN s, GEN Ns)
584 : {
585 28 : x = shallowcopy(x);
586 28 : gel(x,1) = gadd(gel(x,1), s);
587 28 : if (Ns)
588 28 : gel(x,2) = gmul(gel(x,2), Ns);
589 28 : return x;
590 : }
591 :
592 : static GEN
593 14 : poles_translate(GEN x, GEN s, GEN Ns)
594 42 : { pari_APPLY_same(pole_translate(gel(x,i), s, Ns)) }
595 :
596 : /* r / x + O(1) */
597 : static GEN
598 266 : simple_pole(GEN r)
599 : {
600 : GEN S;
601 266 : if (isintzero(r)) return gen_0;
602 217 : S = deg1ser_shallow(gen_0, r, 1);
603 217 : setvalser(S, -1); return S;
604 : }
605 :
606 : GEN
607 98 : lfunshift(GEN ldata, GEN s, long flag, long bitprec)
608 : {
609 98 : pari_sp ltop = avma;
610 : GEN k, k1, L, N, a, b, gam, eps, res;
611 98 : long prec = nbits2prec(bitprec);
612 98 : if (!is_rational_t(typ(s))) pari_err_TYPE("lfunshift",s);
613 98 : ldata = ldata_newprec(lfunmisc_to_ldata_shallow(ldata), prec);
614 98 : a = ldata_get_an(ldata);
615 98 : b = ldata_get_dual(ldata);
616 98 : gam = RgV_Rg_translate(ldata_get_gammavec(ldata), gneg(s));
617 98 : k = gadd(ldata_get_k(ldata), gmul2n(s, 1));
618 98 : k1 = gadd(ldata_get_k1(ldata), s);
619 98 : N = ldata_get_conductor(ldata);
620 98 : eps = ldata_get_rootno(ldata);
621 98 : res = ldata_get_residue(ldata);
622 98 : a = tag(mkvec2(a, s), t_LFUN_SHIFT);
623 98 : if (typ(b) != t_INT)
624 0 : b = tag(mkvec2(b, s), t_LFUN_SHIFT);
625 98 : if (res)
626 98 : switch(typ(res))
627 : {
628 0 : case t_VEC:
629 0 : res = poles_translate(res, s, NULL);
630 0 : break;
631 14 : case t_COL:
632 14 : res = poles_translate(res, s, gpow(N, gmul2n(s, -1), prec));
633 14 : break;
634 84 : default:
635 84 : res = mkvec(mkvec2(gsub(k, s), simple_pole(res)));
636 : }
637 98 : L = mkvecn(res ? 7: 6, a, b, gam, mkvec2(k, k1), N, eps, res);
638 98 : if (flag) L = lfunmul_k(ldata, L, gsub(k, s), bitprec);
639 98 : return gerepilecopy(ltop, L);
640 : }
641 :
642 : /*****************************************************************/
643 : /* L-series from closure */
644 : /*****************************************************************/
645 : static GEN
646 58954 : localfactor(void *E, GEN p, long n)
647 : {
648 58954 : GEN s = closure_callgen2((GEN)E, p, utoi(n));
649 58954 : return direuler_factor(s, n);
650 : }
651 : static GEN
652 1932 : vecan_closure(GEN a, long L, long prec)
653 : {
654 1932 : long ta = typ(a);
655 1932 : GEN gL, Sbad = NULL;
656 :
657 1932 : if (!L) return cgetg(1,t_VEC);
658 1932 : if (ta == t_VEC)
659 : {
660 1015 : long l = lg(a);
661 1015 : if (l == 1) pari_err_TYPE("vecan_closure", a);
662 1015 : ta = typ(gel(a,1));
663 : /* regular vector, return it */
664 1015 : if (ta != t_CLOSURE) return vecslice(a, 1, minss(L,l-1));
665 119 : if (l != 3) pari_err_TYPE("vecan_closure", a);
666 119 : Sbad = gel(a,2);
667 119 : if (typ(Sbad) != t_VEC) pari_err_TYPE("vecan_closure", a);
668 112 : a = gel(a,1);
669 : }
670 917 : else if (ta != t_CLOSURE) pari_err_TYPE("vecan_closure", a);
671 1015 : push_localprec(prec);
672 1015 : gL = stoi(L);
673 1015 : switch(closure_arity(a))
674 : {
675 371 : case 2:
676 371 : a = direuler_bad((void*)a, localfactor, gen_2, gL,gL, Sbad);
677 336 : break;
678 637 : case 1:
679 637 : a = closure_callgen1(a, gL);
680 637 : if (typ(a) != t_VEC) pari_err_TYPE("vecan_closure", a);
681 630 : break;
682 7 : default: pari_err_TYPE("vecan_closure [wrong arity]", a);
683 : a = NULL; /*LCOV_EXCL_LINE*/
684 : }
685 966 : pop_localprec(); return a;
686 : }
687 :
688 : static GEN
689 70 : eulerf_closure(GEN a, GEN p, long prec)
690 : {
691 70 : long ta = typ(a);
692 70 : GEN Sbad = NULL, f;
693 :
694 70 : if (ta == t_VEC)
695 : {
696 0 : long l = lg(a);
697 0 : if (l == 1) pari_err_TYPE("vecan_closure", a);
698 0 : ta = typ(gel(a,1));
699 : /* regular vector, return it */
700 0 : if (ta != t_CLOSURE) return NULL;
701 0 : if (l != 3) pari_err_TYPE("vecan_closure", a);
702 0 : Sbad = gel(a,2);
703 0 : if (typ(Sbad) != t_VEC) pari_err_TYPE("vecan_closure", a);
704 0 : a = gel(a,1);
705 : }
706 70 : else if (ta != t_CLOSURE) pari_err_TYPE("vecan_closure", a);
707 70 : push_localprec(prec);
708 70 : switch(closure_arity(a))
709 : {
710 14 : case 2:
711 14 : f = closure_callgen2(a, p, mkoo()); break;
712 56 : case 1:
713 56 : f = NULL; break;
714 0 : default:
715 0 : f = NULL; pari_err_TYPE("vecan_closure", a);
716 : }
717 70 : pop_localprec(); return f;
718 : }
719 :
720 : /*****************************************************************/
721 : /* L-series of Dirichlet characters. */
722 : /*****************************************************************/
723 :
724 : static GEN
725 4074 : lfunzeta(void)
726 : {
727 4074 : GEN zet = mkvecn(7, NULL, gen_0, NULL, gen_1, gen_1, gen_1, gen_1);
728 4074 : gel(zet,1) = tag(gen_1, t_LFUN_ZETA);
729 4074 : gel(zet,3) = mkvec(gen_0);
730 4074 : return zet;
731 : }
732 :
733 : static GEN
734 3654 : vecan_Kronecker(GEN D, long n)
735 : {
736 3654 : GEN v = cgetg(n+1, t_VECSMALL);
737 3654 : ulong Du = itou_or_0(D);
738 3654 : long i, id, d = Du ? minuu(Du, n): n;
739 42483 : for (i = 1; i <= d; i++) v[i] = krois(D,i);
740 553630 : for (id = i; i <= n; i++,id++) /* periodic mod d */
741 : {
742 549976 : if (id > d) id = 1;
743 549976 : gel(v, i) = gel(v, id);
744 : }
745 3654 : return v;
746 : }
747 :
748 : static GEN
749 6881 : lfunchiquad(GEN D)
750 : {
751 : GEN r;
752 6881 : D = coredisc(D);
753 6881 : if (equali1(D)) return lfunzeta();
754 6055 : if (!isfundamental(D)) pari_err_TYPE("lfunchiquad [not primitive]", D);
755 6055 : r = mkvecn(6, NULL, gen_0, NULL, gen_1, NULL, gen_1);
756 6055 : gel(r,1) = tag(icopy(D), t_LFUN_KRONECKER);
757 6055 : gel(r,3) = mkvec(signe(D) < 0? gen_1: gen_0);
758 6055 : gel(r,5) = mpabs(D);
759 6055 : return r;
760 : }
761 :
762 : /* Begin Hecke characters. Here a character is assumed to be given by a
763 : vector on the generators of the ray class group clgp of CL_m(K).
764 : If clgp = [h,[d1,...,dk],[g1,...,gk]] with dk|...|d2|d1, a character chi
765 : is given by [a1,a2,...,ak] such that chi(gi)=\zeta_di^ai. */
766 :
767 : /* Value of CHI on x, coprime to bnr.mod */
768 : static GEN
769 97293 : chigeneval_i(GEN logx, GEN d, GEN nchi, GEN z, long prec)
770 : {
771 97293 : pari_sp av = avma;
772 97293 : GEN e = FpV_dotproduct(nchi, logx, d);
773 97293 : if (!is_vec_t(typ(z)))
774 1225 : return gerepileupto(av, gpow(z, e, prec));
775 : else
776 : {
777 96068 : ulong i = itou(e);
778 96068 : set_avma(av); return gel(z, i+1);
779 : }
780 : }
781 :
782 : static GEN
783 85050 : chigenevalvec(GEN logx, GEN nchi, GEN z, long prec, long multi)
784 : {
785 85050 : GEN d = gel(nchi,1), x = gel(nchi, 2);
786 85050 : if (multi)
787 35511 : pari_APPLY_same(chigeneval_i(logx, d, gel(x,i), z, prec))
788 : else
789 73416 : return chigeneval_i(logx, d, x, z, prec);
790 : }
791 :
792 : /* return x + yz; y != 0; z = 0,1 "often"; x = 0 "often" */
793 : static GEN
794 1887018 : gaddmul(GEN x, GEN y, GEN z)
795 : {
796 : pari_sp av;
797 1887018 : if (typ(z) == t_INT)
798 : {
799 1684025 : if (!signe(z)) return x;
800 28350 : if (equali1(z)) return gadd(x,y);
801 : }
802 222943 : if (isintzero(x)) return gmul(y,z);
803 119931 : av = avma;
804 119931 : return gerepileupto(av, gadd(x, gmul(y,z)));
805 : }
806 :
807 : static GEN
808 1805916 : gaddmulvec(GEN x, GEN y, GEN z, long multi)
809 : {
810 1805916 : if (multi)
811 244916 : pari_APPLY_same(gaddmul(gel(x,i),gel(y,i),gel(z,i)))
812 : else
813 1724009 : return gaddmul(x,y,z);
814 : }
815 :
816 : static GEN
817 4991 : mkvchi(GEN chi, long n)
818 : {
819 : GEN v;
820 4991 : if (lg(chi) > 1 && is_vec_t(typ(gel(chi,1))))
821 483 : {
822 483 : long d = lg(chi)-1;
823 483 : v = const_vec(n, zerovec(d));
824 483 : gel(v,1) = const_vec(d, gen_1);
825 : }
826 : else
827 4508 : v = vec_ei(n, 1);
828 4991 : return v;
829 : }
830 :
831 : static GEN
832 3647 : vecan_chiZ(GEN an, long n, long prec)
833 : {
834 : forprime_t iter;
835 3647 : GEN G = gel(an,1);
836 3647 : GEN nchi = gel(an,2), gord = gel(nchi,1), chi = gel(nchi,2), z;
837 3647 : GEN gp = cgetipos(3), v = mkvchi(chi, n);
838 3647 : GEN N = znstar_get_N(G);
839 3647 : long ord = itos_or_0(gord);
840 3647 : ulong Nu = itou_or_0(N);
841 3647 : long i, id, d = Nu ? minuu(Nu, n): n;
842 3647 : long multichi= (lg(chi) > 1 && is_vec_t(typ(gel(chi,1))));
843 : ulong p;
844 3647 : if (!multichi && ord && n > (ord>>4))
845 3521 : {
846 3521 : GEN w = ncharvecexpo(G, nchi);
847 3521 : z = grootsof1(ord, prec);
848 37926 : for (i = 1; i <= d; i++)
849 34405 : if (w[i] >= 0) gel(v, i) = gel(z, w[i]+1);
850 : }
851 : else
852 : {
853 126 : z = rootsof1_cx(gord, prec);
854 126 : u_forprime_init(&iter, 2, d);
855 805 : while ((p = u_forprime_next(&iter)))
856 : {
857 : GEN ch;
858 : ulong k;
859 679 : if (!umodiu(N,p)) continue;
860 560 : gp[2] = p;
861 560 : ch = chigenevalvec(znconreylog(G, gp), nchi, z, prec, multichi);
862 560 : gel(v, p) = ch;
863 1582 : for (k = 2*p; k <= (ulong)d; k += p)
864 1022 : gel(v, k) = gaddmulvec(gel(v, k), ch, gel(v, k/p), multichi);
865 : }
866 : }
867 898765 : for (id = i = d+1; i <= n; i++,id++) /* periodic mod d */
868 : {
869 895118 : if (id > d) id = 1;
870 895118 : gel(v, i) = gel(v, id);
871 : }
872 3647 : return v;
873 : }
874 :
875 : static GEN
876 42 : eulerf_chiZ(GEN an, GEN p, long prec)
877 : {
878 42 : GEN G = gel(an,1);
879 42 : GEN nchi = gel(an,2), gord = gel(nchi,1), chi = gel(nchi,2);
880 42 : long multichi= (lg(chi) > 1 && is_vec_t(typ(gel(chi,1))));
881 42 : GEN z = rootsof1_cx(gord, prec);
882 42 : GEN N = znstar_get_N(G);
883 42 : GEN ch = dvdii(N,p) ? gen_0: chigenevalvec(znconreylog(G, p), nchi, z, prec, multichi);
884 42 : return mkrfrac(gen_1, deg1pol_shallow(gneg(ch), gen_1,0));
885 : }
886 :
887 : static GEN
888 1344 : vecan_chigen(GEN an, long n, long prec)
889 : {
890 : forprime_t iter;
891 1344 : GEN bnr = gel(an,1), nf = bnr_get_nf(bnr);
892 1344 : GEN nchi = gel(an,2), gord = gel(nchi,1), chi = gel(nchi,2), z;
893 1344 : GEN gp = cgetipos(3), v = mkvchi(chi, n);
894 1344 : GEN N = gel(bnr_get_mod(bnr), 1), NZ = gcoeff(N,1,1);
895 1344 : long ord = itos_or_0(gord);
896 1344 : long multichi= (lg(chi) > 1 && is_vec_t(typ(gel(chi,1))));
897 : ulong p;
898 :
899 1344 : if (ord && n > (ord>>4))
900 1344 : z = grootsof1(ord, prec);
901 : else
902 0 : z = rootsof1_cx(gord, prec);
903 :
904 1344 : if (nf_get_degree(nf) == 1)
905 : {
906 945 : ulong Nu = itou_or_0(NZ);
907 945 : long i, id, d = Nu ? minuu(Nu, n): n;
908 945 : u_forprime_init(&iter, 2, d);
909 4711 : while ((p = u_forprime_next(&iter)))
910 : {
911 : GEN ch;
912 : ulong k;
913 3766 : if (!umodiu(NZ,p)) continue;
914 2800 : gp[2] = p;
915 2800 : ch = chigenevalvec(isprincipalray(bnr,gp), nchi, z, prec, multichi);
916 2800 : gel(v, p) = ch;
917 7987 : for (k = 2*p; k <= (ulong)d; k += p)
918 5187 : gel(v, k) = gaddmulvec(gel(v, k), ch, gel(v, k/p), multichi);
919 : }
920 11543 : for (id = i = d+1; i <= n; i++,id++) /* periodic mod d */
921 : {
922 10598 : if (id > d) id = 1;
923 10598 : gel(v, i) = gel(v, id);
924 : }
925 : }
926 : else
927 : {
928 399 : GEN BOUND = stoi(n);
929 399 : u_forprime_init(&iter, 2, n);
930 82369 : while ((p = u_forprime_next(&iter)))
931 : {
932 : GEN L;
933 : long j;
934 81970 : int check = !umodiu(NZ,p);
935 81970 : gp[2] = p;
936 81970 : L = idealprimedec_limit_norm(nf, gp, BOUND);
937 163828 : for (j = 1; j < lg(L); j++)
938 : {
939 81858 : GEN pr = gel(L, j), ch;
940 : ulong k, q;
941 81858 : if (check && idealval(nf, N, pr)) continue;
942 81613 : ch = chigenevalvec(isprincipalray(bnr,pr), nchi, z, prec, multichi);
943 81613 : q = upr_norm(pr);
944 81613 : gel(v, q) = gadd(gel(v, q), ch);
945 1881320 : for (k = 2*q; k <= (ulong)n; k += q)
946 1799707 : gel(v, k) = gaddmulvec(gel(v, k), ch, gel(v, k/q), multichi);
947 : }
948 : }
949 : }
950 1344 : return v;
951 : }
952 :
953 : static GEN
954 28 : eulerf_chigen(GEN an, GEN p, long prec)
955 : {
956 28 : GEN bnr = gel(an,1), nf = bnr_get_nf(bnr);
957 28 : GEN nchi = gel(an,2), gord = gel(nchi,1), chi = gel(nchi,2), z;
958 28 : GEN N = gel(bnr_get_mod(bnr), 1), NZ = gcoeff(N,1,1), f;
959 28 : long multichi= (lg(chi) > 1 && is_vec_t(typ(gel(chi,1))));
960 :
961 28 : z = rootsof1_cx(gord, prec);
962 28 : if (nf_get_degree(nf) == 1)
963 : {
964 : GEN ch;
965 0 : if (dvdii(NZ,p)) ch = gen_0;
966 : else
967 : {
968 0 : ch = chigenevalvec(isprincipalray(bnr,p), nchi, z, prec, multichi);
969 0 : if (typ(ch)==t_VEC) return NULL;
970 : }
971 0 : f = deg1pol_shallow(gneg(ch), gen_1, 0);
972 : }
973 : else
974 : {
975 28 : int check = dvdii(NZ,p);
976 28 : GEN L = idealprimedec(nf, p);
977 28 : long j, lL = lg(L);
978 28 : f = pol_1(0);
979 49 : for (j = 1; j < lL; j++)
980 : {
981 35 : GEN pr = gel(L, j), ch;
982 35 : if (check && idealval(nf, N, pr)) ch = gen_0;
983 : else
984 35 : ch = chigenevalvec(isprincipalray(bnr,pr), nchi, z, prec, multichi);
985 35 : if (typ(ch)==t_VEC) return NULL;
986 21 : f = gmul(f, gsub(gen_1, monomial(ch, pr_get_f(pr), 0)));
987 : }
988 : }
989 14 : return mkrfrac(gen_1,f);
990 : }
991 :
992 : static GEN
993 4480 : vec01(long r1, long r2)
994 : {
995 4480 : long d = r1+r2, i;
996 4480 : GEN v = cgetg(d+1,t_VEC);
997 12649 : for (i = 1; i <= r1; i++) gel(v,i) = gen_0;
998 7203 : for ( ; i <= d; i++) gel(v,i) = gen_1;
999 4480 : return v;
1000 : }
1001 :
1002 : /* true nf or t_POL */
1003 : static GEN
1004 1246 : lfunzetak_i(GEN T)
1005 : {
1006 : GEN Vga, N;
1007 : long r1, r2;
1008 1246 : if (typ(T) == t_POL)
1009 : {
1010 651 : T = nfinit0(T, nf_NOLLL, DEFAULTPREC);
1011 651 : if (lg(T) == 3) T = gel(T,1); /* [nf,change of var] */
1012 : }
1013 1246 : if (nf_get_degree(T) == 1) return lfunzeta();
1014 1246 : nf_get_sign(T,&r1,&r2); Vga = vec01(r1+r2,r2);
1015 1246 : N = absi_shallow(nf_get_disc(T));
1016 1246 : return mkvecn(7, tag(T,t_LFUN_NF), gen_0, Vga, gen_1, N, gen_1, gen_0);
1017 : }
1018 : /* truen nf or t_POL */
1019 : static GEN
1020 819 : lfunzetak(GEN T)
1021 819 : { pari_sp av = avma; return gerepilecopy(av, lfunzetak_i(T)); }
1022 :
1023 : /* C = normalized character of order dividing o; renormalize so that it has
1024 : * apparent order o */
1025 : static GEN
1026 609 : char_renormalize(GEN C, GEN o)
1027 : {
1028 609 : GEN oc = gel(C,1), c = gel(C,2);
1029 609 : if (!equalii(o, oc)) c = ZC_Z_mul(c, diviiexact(o, oc));
1030 609 : return c;
1031 : }
1032 : static GEN
1033 224 : vecchar_renormalize(GEN x, GEN o)
1034 833 : { pari_APPLY_same(char_renormalize(gel(x,i), o)); }
1035 :
1036 : /* G is a bid of nftyp typ_BIDZ */
1037 : static GEN
1038 8421 : lfunchiZ(GEN G, GEN CHI)
1039 : {
1040 8421 : pari_sp av = avma;
1041 8421 : GEN sig = NULL, N = bid_get_ideal(G), nchi, r;
1042 : int real;
1043 : long s;
1044 :
1045 8421 : if (typ(N) != t_INT) pari_err_TYPE("lfunchiZ", G);
1046 8421 : if (typ(CHI) == t_VEC && !RgV_is_ZV(CHI))
1047 21 : {
1048 42 : GEN C, G0 = G, o = gen_1;
1049 42 : long i, l = lg(CHI);
1050 42 : nchi = cgetg(l, t_VEC);
1051 42 : N = znconreyconductor(G, gel(CHI,1), &C);
1052 35 : if (typ(N) != t_INT) G = znstar0(N, 1);
1053 35 : s = zncharisodd(G, C);
1054 91 : for (i = 1; i < l; i++)
1055 : {
1056 70 : if (i > 1)
1057 : {
1058 35 : if (!gequal(N, znconreyconductor(G0, gel(CHI,i), &C))
1059 28 : || zncharisodd(G, C) != s)
1060 14 : pari_err_TYPE("lfuncreate [different conductors]", CHI);
1061 : }
1062 56 : C = znconreylog_normalize(G, C);
1063 56 : o = lcmii(o, gel(C,1)); /* lcm with charorder */
1064 56 : gel(nchi,i) = C;
1065 : }
1066 21 : nchi = mkvec2(o, vecchar_renormalize(nchi, o));
1067 21 : if (typ(N) != t_INT) N = gel(N,1);
1068 : }
1069 : else
1070 : {
1071 8379 : N = znconreyconductor(G, CHI, &CHI);
1072 8379 : if (typ(N) != t_INT)
1073 : {
1074 7 : if (equali1(gel(N,1))) { set_avma(av); return lfunzeta(); }
1075 0 : G = znstar0(N, 1);
1076 0 : N = gel(N,1);
1077 : }
1078 : /* CHI now primitive on G */
1079 8372 : switch(itou_or_0(zncharorder(G, CHI)))
1080 : {
1081 2415 : case 1: set_avma(av); return lfunzeta();
1082 2982 : case 2: if (zncharisodd(G,CHI)) N = negi(N);
1083 2982 : return gerepileupto(av, lfunchiquad(N));
1084 : }
1085 2975 : nchi = znconreylog_normalize(G, CHI);
1086 2975 : s = zncharisodd(G, CHI);
1087 : }
1088 2996 : sig = mkvec(s? gen_1: gen_0);
1089 2996 : real = abscmpiu(gel(nchi,1), 2) <= 0;
1090 2996 : r = mkvecn(6, tag(mkvec2(G,nchi), t_LFUN_CHIZ),
1091 : real? gen_0: gen_1, sig, gen_1, N, gen_0);
1092 2996 : return gerepilecopy(av, r);
1093 : }
1094 :
1095 : static GEN
1096 1309 : lfunchigen(GEN bnr, GEN CHI)
1097 : {
1098 1309 : pari_sp av = avma;
1099 : GEN N, sig, Ldchi, nf, nchi, NN;
1100 : long r1, r2, n1;
1101 : int real;
1102 :
1103 1309 : if (typ(CHI) == t_VEC && !RgV_is_ZV(CHI))
1104 203 : {
1105 210 : long map, i, l = lg(CHI);
1106 210 : GEN bnr0 = bnr, D, chi = gel(CHI,1), o = gen_1;
1107 210 : nchi = cgetg(l, t_VEC);
1108 210 : bnr_char_sanitize(&bnr, &chi);
1109 210 : D = cyc_normalize(bnr_get_cyc(bnr));
1110 210 : N = bnr_get_mod(bnr);
1111 210 : map = (bnr != bnr0);
1112 784 : for (i = 1; i < l; i++)
1113 : {
1114 581 : if (i > 1)
1115 : {
1116 371 : chi = gel(CHI,i);
1117 371 : if (!map)
1118 : {
1119 350 : if (!bnrisconductor(bnr, chi))
1120 7 : pari_err_TYPE("lfuncreate [different conductors]", CHI);
1121 : }
1122 : else
1123 : {
1124 21 : if (!gequal(bnrconductor_raw(bnr0, chi), N))
1125 0 : pari_err_TYPE("lfuncreate [different conductors]", CHI);
1126 21 : chi = bnrchar_primitive_raw(bnr0, bnr, chi);
1127 : }
1128 : }
1129 574 : chi = char_normalize(chi, D);
1130 574 : o = lcmii(o, gel(chi,1)); /* lcm with charorder */
1131 574 : gel(nchi,i) = chi;
1132 : }
1133 203 : nchi = mkvec2(o, vecchar_renormalize(nchi, o));
1134 : }
1135 : else
1136 : {
1137 1099 : bnr_char_sanitize(&bnr, &CHI);
1138 1099 : nchi = NULL; /* now CHI is primitive wrt bnr */
1139 : }
1140 :
1141 1302 : N = bnr_get_mod(bnr);
1142 1302 : nf = bnr_get_nf(bnr);
1143 1302 : n1 = lg(vec01_to_indices(gel(N,2))) - 1; /* vecsum(N[2]) */
1144 1302 : N = gel(N,1);
1145 1302 : NN = mulii(idealnorm(nf, N), absi_shallow(nf_get_disc(nf)));
1146 1302 : if (!nchi)
1147 : {
1148 1099 : if (equali1(NN)) { set_avma(av); return lfunzeta(); }
1149 658 : if (ZV_equal0(CHI)) return gerepilecopy(av, lfunzetak_i(bnr_get_nf(bnr)));
1150 581 : nchi = char_normalize(CHI, cyc_normalize(bnr_get_cyc(bnr)));
1151 : }
1152 784 : real = abscmpiu(gel(nchi,1), 2) <= 0;
1153 784 : nf_get_sign(nf, &r1, &r2);
1154 784 : sig = vec01(r1+r2-n1, r2+n1);
1155 784 : Ldchi = mkvecn(6, tag(mkvec2(bnr, nchi), t_LFUN_CHIGEN),
1156 : real? gen_0: gen_1, sig, gen_1, NN, gen_0);
1157 784 : return gerepilecopy(av, Ldchi);
1158 : }
1159 :
1160 : /* Find all characters of clgp whose kernel contain group given by HNF H.
1161 : * Set *pcnj[i] to the conductor */
1162 : static GEN
1163 518 : chigenkerfind(GEN bnr, GEN H, GEN *pcnj)
1164 : {
1165 518 : GEN res, cnj, L = bnrchar(bnr, H, NULL);
1166 518 : long i, k, l = lg(L);
1167 :
1168 518 : res = cgetg(l, t_VEC);
1169 518 : *pcnj = cnj = cgetg(l, t_VEC);
1170 1911 : for (i = k = 1; i < l; i++)
1171 : {
1172 1393 : GEN chi = gel(L,i);
1173 1393 : gel(res, k) = chi;
1174 1393 : gel(cnj, k) = ZV_equal0(chi)? gen_0: bnrconductor_raw(bnr, chi);
1175 1393 : k++;
1176 : }
1177 518 : setlg(cnj, k);
1178 518 : setlg(res, k); return res;
1179 : }
1180 :
1181 : static GEN
1182 518 : vec_classes(GEN A, GEN F)
1183 : {
1184 518 : GEN w = vec_equiv(F);
1185 518 : long i, l = lg(w);
1186 518 : GEN V = cgetg(l, t_VEC);
1187 1575 : for (i = 1; i < l; i++) gel(V,i) = vecpermute(A,gel(w,i));
1188 518 : return V;
1189 : }
1190 :
1191 : static GEN
1192 990451 : abelrel_pfactor(GEN bnr, GEN pr, GEN U, GEN D, GEN h)
1193 : {
1194 990451 : GEN v = bnrisprincipalmod(bnr, pr, h, 0);
1195 990451 : GEN E = ZV_ZV_mod(ZM_ZC_mul(U, v), D);
1196 990451 : ulong o = itou(charorder(D, E)), f = pr_get_f(pr);
1197 990451 : return gpowgs(gsub(gen_1, monomial(gen_1, f * o, 0)), itou(h) / o);
1198 : }
1199 :
1200 : static GEN
1201 687925 : abelrel_factor(GEN bnr, GEN C, GEN p, GEN mod, GEN U, GEN D, GEN h)
1202 : {
1203 687925 : GEN nf = bnr_get_nf(bnr), F = pol_1(0), prid = idealprimedec(nf,p);
1204 687925 : GEN mod2 = shallowcopy(mod);
1205 687925 : long i, l = lg(prid);
1206 1678376 : for (i = 1; i < l; i++)
1207 : {
1208 990451 : GEN pr = gel(prid, i), Fpr;
1209 990451 : long v = idealval(nf,mod,pr);
1210 990451 : if (v > 0)
1211 : {
1212 : GEN bnr2, C2, U2, D2, h2;
1213 161 : gel(mod2, 1) = idealdivpowprime(nf, gel(mod, 1), pr, utoi(v));
1214 161 : bnr2 = bnrinitmod(bnr, mod2, 0, h);
1215 161 : C2 = bnrmap(bnrmap(bnr, bnr2), C);
1216 161 : D2 = ZM_snfall_i(C2, &U2, NULL, 1);
1217 161 : h2 = ZV_prod(D2);
1218 161 : Fpr = abelrel_pfactor(bnr2, pr, U2, D2, h2);
1219 : }
1220 : else
1221 990290 : Fpr = abelrel_pfactor(bnr, pr, U, D, h);
1222 990451 : F = ZX_mul(F, Fpr);
1223 : }
1224 687925 : return gcopy(mkrfrac(gen_1, F));
1225 : }
1226 :
1227 : static GEN
1228 42 : eulerf_abelrel(GEN an, GEN p)
1229 : {
1230 42 : GEN bnr = gel(an,1), C = gel(an,2), mod = gel(an,3);
1231 42 : GEN U, D = ZM_snfall_i(C, &U, NULL, 1), h = ZV_prod(D);
1232 42 : return abelrel_factor(bnr, C, p, mod, U, D, h);
1233 : }
1234 :
1235 : struct direuler_abelrel
1236 : {
1237 : GEN bnr, C, mod, U, D, h;
1238 : };
1239 :
1240 : static GEN
1241 687883 : _direuler_abelrel(void *E, GEN p)
1242 : {
1243 687883 : struct direuler_abelrel *s = (struct direuler_abelrel*) E;
1244 687883 : return abelrel_factor(s->bnr, s->C, p, s->mod, s->U, s->D, s->h);
1245 : }
1246 :
1247 : static GEN
1248 168 : vecan_abelrel(GEN an, long N)
1249 : {
1250 : struct direuler_abelrel s;
1251 168 : s.bnr = gel(an,1);
1252 168 : s.C = gel(an,2);
1253 168 : s.mod = gel(an,3);
1254 168 : s.D = ZM_snfall_i(s.C, &s.U, NULL, 1);
1255 168 : s.h = ZV_prod(s.D);
1256 168 : return direuler((void*)&s, _direuler_abelrel, gen_1, stoi(N), NULL);
1257 : }
1258 :
1259 : static GEN
1260 539 : lfunabelrel_i(GEN bnr, GEN H, GEN mod)
1261 : {
1262 539 : GEN NrD = bnrdisc(bnr, H, 0), N = absi_shallow(gel(NrD,3));
1263 539 : long n = itos(gel(NrD,1)), r1 = itos(gel(NrD,2)), r2 = (n-r1)>>1;
1264 539 : if (!mod) mod = bnrconductor(bnr, H, 0);
1265 539 : return mkvecn(7, tag(mkvec3(bnr,H,mod),t_LFUN_ABELREL),
1266 : gen_0, vec01(r1+r2, r2), gen_1, N, gen_1, gen_0);
1267 : }
1268 : static GEN
1269 21 : lfunabelrel(GEN bnr, GEN H, GEN mod)
1270 21 : { pari_sp av = avma; return gerepilecopy(av, lfunabelrel_i(bnr, H, mod)); }
1271 :
1272 :
1273 : static GEN
1274 1057 : lfunchiinit(GEN bnr, GEN chi, GEN dom, long der, long bitprec)
1275 : {
1276 1057 : GEN L = lfunchigen(bnr, lg(chi)==2 ? gel(chi,1): chi);
1277 1057 : return lfuninit(L, dom, der, bitprec);
1278 : }
1279 : static GEN
1280 518 : veclfunchiinit(GEN bnr, GEN x, GEN dom, long der, long bitprec)
1281 1575 : { pari_APPLY_same(lfunchiinit(bnr, gel(x,i), dom, der, bitprec)); }
1282 : GEN
1283 518 : lfunabelianrelinit(GEN bnr, GEN H, GEN dom, long der, long bitprec)
1284 : {
1285 518 : GEN cnj, M, D, C = chigenkerfind(bnr, H, &cnj);
1286 : long l;
1287 518 : C = vec_classes(C, cnj); l = lg(C);
1288 518 : M = mkvec3(veclfunchiinit(bnr, C, dom, der, bitprec),
1289 : const_vecsmall(l-1, 1), const_vecsmall(l-1, 0));
1290 518 : D = mkvec2(dom, mkvecsmall2(der, bitprec));
1291 518 : return lfuninit_make(t_LDESC_PRODUCT, lfunabelrel_i(bnr, H, NULL), M, D);
1292 : }
1293 :
1294 : /*****************************************************************/
1295 : /* Dedekind zeta functions */
1296 : /*****************************************************************/
1297 : /* true nf */
1298 : static GEN
1299 2107 : dirzetak0(GEN nf, ulong N)
1300 : {
1301 2107 : GEN vect, c, c2, T = nf_get_pol(nf), index = nf_get_index(nf);
1302 2107 : pari_sp av = avma, av2;
1303 2107 : const ulong SQRTN = usqrt(N);
1304 : ulong i, p, lx;
1305 2107 : long court[] = {evaltyp(t_INT)|_evallg(3), evalsigne(1)|evallgefint(3),0};
1306 : forprime_t S;
1307 :
1308 2107 : c = cgetalloc(N+1, t_VECSMALL);
1309 2107 : c2 = cgetalloc(N+1, t_VECSMALL);
1310 2671991 : c2[1] = c[1] = 1; for (i=2; i<=N; i++) c[i] = 0;
1311 2107 : u_forprime_init(&S, 2, N); av2 = avma;
1312 348579 : while ( (p = u_forprime_next(&S)) )
1313 : {
1314 346472 : set_avma(av2);
1315 346472 : if (umodiu(index, p)) /* p does not divide index */
1316 346108 : vect = gel(Flx_degfact(ZX_to_Flx(T,p), p),1);
1317 : else
1318 : {
1319 364 : court[2] = p;
1320 364 : vect = idealprimedec_degrees(nf,court);
1321 : }
1322 346472 : lx = lg(vect);
1323 346472 : if (p <= SQRTN)
1324 35245 : for (i=1; i<lx; i++)
1325 : {
1326 23548 : ulong qn, q = upowuu(p, vect[i]); /* Norm P[i] */
1327 23548 : if (!q || q > N) break;
1328 20797 : memcpy(c2 + 2, c + 2, (N-1)*sizeof(long));
1329 : /* c2[i] <- c[i] + sum_{k = 1}^{v_q(i)} c[i/q^k] for all i <= N */
1330 42889 : for (qn = q; qn <= N; qn *= q)
1331 : {
1332 42889 : ulong k0 = N/qn, k, k2; /* k2 = k*qn */
1333 3760050 : for (k = k0, k2 = k*qn; k > 0; k--, k2 -=qn) c2[k2] += c[k];
1334 42889 : if (q > k0) break; /* <=> q*qn > N */
1335 : }
1336 20797 : swap(c, c2);
1337 : }
1338 : else /* p > sqrt(N): simpler */
1339 656292 : for (i=1; i<lx; i++)
1340 : {
1341 : ulong k, k2; /* k2 = k*p */
1342 581798 : if (vect[i] > 1) break;
1343 : /* c2[i] <- c[i] + sum_{k = 1}^{v_q(i)} c[i/q^k] for all i <= N */
1344 1906723 : for (k = N/p, k2 = k*p; k > 0; k--, k2 -= p) c[k2] += c[k];
1345 : }
1346 : }
1347 2107 : set_avma(av);
1348 2107 : pari_free(c2); return c;
1349 : }
1350 :
1351 : static GEN
1352 168 : eulerf_zetak(GEN nf, GEN p)
1353 : {
1354 168 : GEN v, f = pol_1(0);
1355 : long i, l;
1356 168 : if (dvdii(nf_get_index(nf), p)) /* p does not divide index */
1357 7 : v = idealprimedec_degrees(nf,p);
1358 : else
1359 161 : v = gel(FpX_degfact(nf_get_pol(nf), p), 1);
1360 168 : l = lg(v);
1361 406 : for (i = 1; i < l; i++) f = ZX_sub(f, RgX_shift_shallow(f, v[i]));
1362 168 : retmkrfrac(gen_1, ZX_copy(f));
1363 : }
1364 :
1365 : GEN
1366 2107 : dirzetak(GEN nf, GEN b)
1367 : {
1368 : GEN z, c;
1369 : long n;
1370 :
1371 2107 : if (typ(b) != t_INT) pari_err_TYPE("dirzetak",b);
1372 2107 : if (signe(b) <= 0) return cgetg(1,t_VEC);
1373 2107 : nf = checknf(nf);
1374 2107 : n = itou_or_0(b); if (!n) pari_err_OVERFLOW("dirzetak");
1375 2107 : c = dirzetak0(nf, n);
1376 2107 : z = vecsmall_to_vec(c); pari_free(c); return z;
1377 : }
1378 :
1379 : static GEN
1380 658 : linit_get_mat(GEN linit)
1381 : {
1382 658 : if (linit_get_type(linit)==t_LDESC_PRODUCT)
1383 161 : return lfunprod_get_fact(linit_get_tech(linit));
1384 : else
1385 497 : return mkvec3(mkvec(linit), mkvecsmall(1), mkvecsmall(0));
1386 : }
1387 :
1388 : static GEN
1389 329 : lfunproduct(GEN ldata, GEN linit1, GEN linit2, GEN domain)
1390 : {
1391 329 : GEN M1 = linit_get_mat(linit1);
1392 329 : GEN M2 = linit_get_mat(linit2);
1393 329 : GEN M3 = mkvec3(shallowconcat(gel(M1, 1), gel(M2, 1)),
1394 329 : vecsmall_concat(gel(M1, 2), gel(M2, 2)),
1395 329 : vecsmall_concat(gel(M1, 3), gel(M2, 3)));
1396 329 : return lfuninit_make(t_LDESC_PRODUCT, ldata, M3, domain);
1397 : }
1398 : static GEN lfunzetakinit_artin(GEN nf, GEN gal, GEN dom, long der, long bit);
1399 : /* true nf */
1400 : static GEN
1401 329 : lfunzetakinit_quotient(GEN nf, GEN polk, GEN dom, long der, long bitprec)
1402 : {
1403 : GEN ak, an, nfk, Vga, ldata, N, Lk, LKk, domain;
1404 : long r1k, r2k, r1, r2;
1405 :
1406 329 : nf_get_sign(nf,&r1,&r2);
1407 329 : nfk = nfinit(polk, nbits2prec(bitprec));
1408 329 : Lk = lfunzetakinit(nfk, dom, der, bitprec); /* zeta_k */
1409 329 : nf_get_sign(nfk,&r1k,&r2k);
1410 329 : Vga = vec01((r1+r2) - (r1k+r2k), r2-r2k);
1411 329 : N = absi_shallow(diviiexact(nf_get_disc(nf), nf_get_disc(nfk)));
1412 329 : ak = nf_get_degree(nf)==1 ? tag(gen_1, t_LFUN_ZETA): tag(nfk, t_LFUN_NF);
1413 329 : an = tag(mkvec2(tag(nf,t_LFUN_NF), ak), t_LFUN_DIV);
1414 329 : ldata = mkvecn(6, an, gen_0, Vga, gen_1, N, gen_1);
1415 329 : LKk = lfuninit(ldata, dom, der, bitprec); /* zeta_K/zeta_k */
1416 329 : domain = mkvec2(dom, mkvecsmall2(der, bitprec));
1417 329 : return lfunproduct(lfunzetak_i(nf), Lk, LKk, domain);
1418 : }
1419 : /* true nf */
1420 : GEN
1421 945 : lfunzetakinit(GEN nf, GEN dom, long der, long bitprec)
1422 : {
1423 945 : long n, d = nf_get_degree(nf);
1424 945 : GEN L, Q, R, T = nf_get_pol(nf);
1425 945 : if (d == 1) return lfuninit(lfunzeta(), dom, der, bitprec);
1426 777 : if (d > 2)
1427 : {
1428 441 : GEN G = galoisinit(nf, NULL);
1429 441 : if (isintzero(G))
1430 : { /* not Galois */
1431 329 : GEN S = nfsubfields(nf, 0); n = lg(S)-1;
1432 329 : return lfunzetakinit_quotient(nf, gmael(S,n-1,1), dom, der, bitprec);
1433 : }
1434 112 : if (!group_isabelian(galois_group(G))) /* Galois, not Abelian */
1435 21 : return lfunzetakinit_artin(nf, G, dom, der, bitprec);
1436 : }
1437 : /* Abelian over Q */
1438 427 : Q = Buchall(pol_x(1), 0, nbits2prec(bitprec));
1439 427 : T = shallowcopy(T); setvarn(T,0);
1440 427 : R = rnfconductor0(Q, T, 1);
1441 427 : L = lfunabelianrelinit(gel(R,2), gel(R,3), dom, der, bitprec);
1442 427 : delete_var(); return L;
1443 : }
1444 :
1445 : /***************************************************************/
1446 : /* Elliptic Curves and Modular Forms */
1447 : /***************************************************************/
1448 :
1449 : static GEN
1450 203 : lfunellnf(GEN e)
1451 : {
1452 203 : pari_sp av = avma;
1453 203 : GEN ldata = cgetg(7, t_VEC), nf = ellnf_get_nf(e);
1454 203 : GEN N = gel(ellglobalred(e), 1);
1455 203 : long n = nf_get_degree(nf);
1456 203 : gel(ldata, 1) = tag(e, t_LFUN_ELL);
1457 203 : gel(ldata, 2) = gen_0;
1458 203 : gel(ldata, 3) = vec01(n, n);
1459 203 : gel(ldata, 4) = gen_2;
1460 203 : gel(ldata, 5) = mulii(idealnorm(nf,N), sqri(nf_get_disc(nf)));
1461 203 : gel(ldata, 6) = stoi(ellrootno_global(e));
1462 203 : return gerepilecopy(av, ldata);
1463 : }
1464 :
1465 : static GEN
1466 3850 : lfunellQ(GEN e)
1467 : {
1468 3850 : pari_sp av = avma;
1469 3850 : GEN ldata = cgetg(7, t_VEC);
1470 3850 : gel(ldata, 1) = tag(ellanal_globalred(e, NULL), t_LFUN_ELL);
1471 3850 : gel(ldata, 2) = gen_0;
1472 3850 : gel(ldata, 3) = mkvec2(gen_0, gen_1);
1473 3850 : gel(ldata, 4) = gen_2;
1474 3850 : gel(ldata, 5) = ellQ_get_N(e);
1475 3850 : gel(ldata, 6) = stoi(ellrootno_global(e));
1476 3850 : return gerepilecopy(av, ldata); /* ellanal_globalred not gerepile-safe */
1477 : }
1478 :
1479 : static GEN
1480 4053 : lfunell(GEN e)
1481 : {
1482 4053 : long t = ell_get_type(e);
1483 4053 : switch(t)
1484 : {
1485 3850 : case t_ELL_Q: return lfunellQ(e);
1486 203 : case t_ELL_NF:return lfunellnf(e);
1487 : }
1488 0 : pari_err_TYPE("lfun",e);
1489 : return NULL; /*LCOV_EXCL_LINE*/
1490 : }
1491 :
1492 : static GEN
1493 140 : ellsympow_gamma(long m)
1494 : {
1495 140 : GEN V = cgetg(m+2, t_VEC);
1496 140 : long i = 1, j;
1497 140 : if (!odd(m)) gel(V, i++) = stoi(-2*(m>>2));
1498 364 : for (j = (m+1)>>1; j > 0; i+=2, j--)
1499 : {
1500 224 : gel(V,i) = stoi(1-j);
1501 224 : gel(V,i+1) = stoi(1-j+1);
1502 : }
1503 140 : return V;
1504 : }
1505 :
1506 : static GEN
1507 86591 : ellsympow_trace(GEN p, GEN t, long m)
1508 : {
1509 86591 : long k, n = m >> 1;
1510 86591 : GEN tp = gpowers0(sqri(t), n, odd(m)? t: NULL);
1511 86866 : GEN pp = gen_1, b = gen_1, r = gel(tp,n+1);
1512 245264 : for(k=1; k<=n; k++)
1513 : {
1514 : GEN s;
1515 158853 : pp = mulii(pp, p);
1516 158510 : b = diviuexact(muliu(b, (m-(2*k-1))*(m-(2*k-2))), k*(m-(k-1)));
1517 157737 : s = mulii(mulii(b, gel(tp,1+n-k)), pp);
1518 158152 : r = odd(k) ? subii(r, s): addii(r, s);
1519 : }
1520 86411 : return r;
1521 : }
1522 :
1523 : static GEN
1524 3129 : ellsympow_abelian(GEN p, GEN ap, long m, long o)
1525 : {
1526 3129 : pari_sp av = avma;
1527 3129 : long i, M, n = (m+1)>>1;
1528 : GEN pk, tv, pn, pm, F, v;
1529 3129 : if (!odd(o))
1530 : {
1531 0 : if (odd(m)) return pol_1(0);
1532 0 : M = m >> 1; o >>= 1;
1533 : }
1534 : else
1535 3129 : M = m * ((o+1) >> 1);
1536 3129 : pk = gpowers(p,n); pn = gel(pk,n+1);
1537 3129 : tv = cgetg(m+2,t_VEC);
1538 3129 : gel(tv, 1) = gen_2;
1539 3129 : gel(tv, 2) = ap;
1540 10668 : for (i = 3; i <= m+1; i++)
1541 7539 : gel(tv,i) = subii(mulii(ap,gel(tv,i-1)), mulii(p,gel(tv,i-2)));
1542 3129 : pm = odd(m)? mulii(gel(pk,n), pn): sqri(pn); /* cheap p^m */
1543 3129 : F = deg2pol_shallow(pm, gen_0, gen_1, 0);
1544 3129 : v = odd(m) ? pol_1(0): deg1pol_shallow(negi(pn), gen_1, 0);
1545 9450 : for (i = M % o; i < n; i += o) /* o | m-2*i */
1546 : {
1547 6321 : gel(F,3) = negi(mulii(gel(tv,m-2*i+1), gel(pk,i+1)));
1548 6321 : v = ZX_mul(v, F);
1549 : }
1550 3129 : return gerepilecopy(av, v);
1551 : }
1552 :
1553 : static GEN
1554 89811 : ellsympow(GEN E, ulong m, GEN p, long n)
1555 : {
1556 89811 : pari_sp av = avma;
1557 89811 : GEN ap = ellap(E, p);
1558 89633 : if (n <= 2)
1559 : {
1560 86534 : GEN t = ellsympow_trace(p, ap, m);
1561 86414 : return deg1pol_shallow(t, gen_1, 0);
1562 : }
1563 : else
1564 3099 : return gerepileupto(av, RgXn_inv_i(ellsympow_abelian(p, ap, m, 1), n));
1565 : }
1566 :
1567 : GEN
1568 5656 : direllsympow_worker(GEN P, ulong X, GEN E, ulong m)
1569 : {
1570 5656 : pari_sp av = avma;
1571 5656 : long i, l = lg(P);
1572 5656 : GEN W = cgetg(l, t_VEC);
1573 95454 : for(i = 1; i < l; i++)
1574 : {
1575 89798 : ulong p = uel(P,i);
1576 89798 : long d = ulogint(X, p) + 1; /* minimal d such that p^d > X */
1577 89809 : gel(W,i) = ellsympow(E, m, utoi(uel(P,i)), d);
1578 : }
1579 5656 : return gerepilecopy(av, mkvec2(P,W));
1580 : }
1581 :
1582 : static GEN
1583 56 : eulerf_bad(GEN bad, GEN p)
1584 : {
1585 56 : long i, l = lg(bad);
1586 112 : for (i = 1; i < l; i++)
1587 56 : if (equalii(gmael(bad,i,1), p))
1588 0 : return gmael(bad,i,2);
1589 56 : return NULL;
1590 : }
1591 :
1592 : static GEN
1593 343 : vecan_ellsympow(GEN an, long n)
1594 : {
1595 343 : GEN nn = utoi(n), crvm = gel(an,1), bad = gel(an,2);
1596 343 : GEN worker = snm_closure(is_entry("_direllsympow_worker"), crvm);
1597 343 : return pardireuler(worker, gen_2, nn, nn, bad);
1598 : }
1599 :
1600 : static GEN
1601 28 : eulerf_ellsympow(GEN an, GEN p)
1602 : {
1603 28 : GEN crvm = gel(an,1), bad = gel(an,2), E = gel(crvm,1);
1604 28 : GEN f = eulerf_bad(bad, p);
1605 28 : if (f) return f;
1606 28 : retmkrfrac(gen_1,ellsympow_abelian(p, ellap(E, p), itos(gel(crvm,2)), 1));
1607 : }
1608 :
1609 : static long
1610 196 : ellsympow_betam(long o, long m)
1611 : {
1612 196 : const long c3[]={3, -1, 1};
1613 196 : const long c12[]={6, -2, 2, 0, 4, -4};
1614 196 : const long c24[]={12, -2, -4, 6, 4, -10};
1615 196 : if (!odd(o) && odd(m)) return 0;
1616 161 : switch(o)
1617 : {
1618 0 : case 1: return m+1;
1619 14 : case 2: return m+1;
1620 84 : case 3: case 6: return (m+c3[m%3])/3;
1621 0 : case 4: return m%4 == 0 ? (m+2)/2: m/2;
1622 21 : case 8: return m%4 == 0 ? (m+4)/4: (m-2)/4;
1623 35 : case 12: return (m+c12[(m%12)/2])/6;
1624 7 : case 24: return (m+c24[(m%12)/2])/12;
1625 : }
1626 0 : return 0;
1627 : }
1628 :
1629 : static long
1630 98 : ellsympow_epsm(long o, long m) { return m + 1 - ellsympow_betam(o, m); }
1631 :
1632 : static GEN
1633 98 : ellsympow_multred(GEN E, GEN p, long m, long vN, long *cnd, long *w)
1634 : {
1635 98 : if (vN == 1 || !odd(m))
1636 : {
1637 98 : GEN s = (odd(m) && signe(ellap(E,p)) < 0)? gen_1: gen_m1;
1638 98 : *cnd = m;
1639 98 : *w = odd(m)? ellrootno(E, p): 1;
1640 98 : return deg1pol_shallow(s, gen_1, 0);
1641 : }
1642 : else
1643 : {
1644 0 : *cnd = equaliu(p,2)? ((m+1)>>1) * vN: m+1;
1645 0 : *w = (m & 3) == 1? ellrootno(E, p): 1;
1646 0 : return pol_1(0);
1647 : }
1648 : }
1649 :
1650 : static GEN
1651 98 : ellsympow_nonabelian(GEN p, long m, long bet)
1652 : {
1653 98 : GEN q = powiu(p, m >> 1), q2 = sqri(q), F;
1654 98 : if (odd(m))
1655 : {
1656 35 : q2 = mulii(q2, p); /* p^m */
1657 35 : return gpowgs(deg2pol_shallow(q2, gen_0, gen_1, 0), bet>>1);
1658 : }
1659 63 : togglesign_safe(&q2);
1660 63 : F = gpowgs(deg2pol_shallow(q2, gen_0, gen_1, 0), bet>>1);
1661 63 : if (!odd(bet)) return F;
1662 28 : if (m%4 != 2) togglesign_safe(&q);
1663 28 : return gmul(F, deg1pol_shallow(q, gen_1, 0));
1664 : }
1665 :
1666 : static long
1667 0 : safe_Z_pvalrem(GEN n, GEN p, GEN *pr)
1668 0 : { return signe(n)==0? -1: Z_pvalrem(n, p, pr); }
1669 :
1670 : static GEN
1671 0 : c4c6_ap(GEN c4, GEN c6, GEN p)
1672 : {
1673 0 : GEN N = Fp_ellcard(Fp_muls(c4, -27, p), Fp_muls(c6, -54, p), p);
1674 0 : return subii(addiu(p, 1), N);
1675 : }
1676 :
1677 : static GEN
1678 0 : ellsympow_abelian_twist(GEN E, GEN p, long m, long o)
1679 : {
1680 0 : GEN ap, c4t, c6t, c4 = ell_get_c4(E), c6 = ell_get_c6(E);
1681 0 : long v4 = safe_Z_pvalrem(c4, p, &c4t);
1682 0 : long v6 = safe_Z_pvalrem(c6, p, &c6t);
1683 0 : if (v6>=0 && (v4==-1 || 3*v4>=2*v6)) c6 = c6t;
1684 0 : if (v4>=0 && (v6==-1 || 3*v4<=2*v6)) c4 = c4t;
1685 0 : ap = c4c6_ap(c4, c6, p);
1686 0 : return ellsympow_abelian(p, ap, m, o);
1687 : }
1688 :
1689 : static GEN
1690 0 : ellsympow_goodred(GEN E, GEN p, long m, long *cnd, long *w)
1691 : {
1692 0 : long o = 12/cgcd(12, Z_pval(ell_get_disc(E), p));
1693 0 : long bet = ellsympow_betam(o, m);
1694 0 : long eps = m + 1 - bet;
1695 0 : *w = odd(m) && odd(eps>>1) ? ellrootno(E,p): 1;
1696 0 : *cnd = eps;
1697 0 : if (umodiu(p, o) == 1)
1698 0 : return ellsympow_abelian_twist(E, p, m, o);
1699 : else
1700 0 : return ellsympow_nonabelian(p, m, bet);
1701 : }
1702 :
1703 : static long
1704 70 : ellsympow_inertia3(GEN E, long vN)
1705 : {
1706 70 : long vD = Z_lval(ell_get_disc(E), 3);
1707 70 : if (vN==2) return vD%2==0 ? 2: 4;
1708 70 : if (vN==4) return vD%4==0 ? 3: 6;
1709 70 : if (vN==3 || vN==5) return 12;
1710 0 : return 0;
1711 : }
1712 :
1713 : static long
1714 70 : ellsympow_deltam3(long o, long m, long vN)
1715 : {
1716 70 : if (o==3 || o==6) return ellsympow_epsm(3, m);
1717 70 : if (o==12 && vN ==3) return (ellsympow_epsm(3, m))/2;
1718 0 : if (o==12 && vN ==5) return (ellsympow_epsm(3, m))*3/2;
1719 0 : return 0;
1720 : }
1721 :
1722 : static long
1723 0 : ellsympow_isabelian3(GEN E)
1724 : {
1725 0 : ulong c4 = umodiu(ell_get_c4(E),81), c6 = umodiu(ell_get_c6(E), 243);
1726 0 : return (c4 == 27 || (c4%27==9 && (c6==108 || c6==135)));
1727 : }
1728 :
1729 : static long
1730 35 : ellsympow_rootno3(GEN E, GEN p, long o, long m)
1731 : {
1732 35 : const long w6p[]={1,-1,-1,-1,1,1};
1733 35 : const long w6n[]={-1,1,-1,1,-1,1};
1734 35 : const long w12p[]={1,1,-1,1,1,1};
1735 35 : const long w12n[]={-1,-1,-1,-1,-1,1};
1736 35 : long w = ellrootno(E, p), mm = (m%12)>>1;
1737 35 : switch(o)
1738 : {
1739 0 : case 2: return m%4== 1 ? -1: 1;
1740 0 : case 6: return w == 1 ? w6p[mm]: w6n[mm];
1741 35 : case 12: return w == 1 ? w12p[mm]: w12n[mm];
1742 0 : default: return 1;
1743 : }
1744 : }
1745 :
1746 : static GEN
1747 70 : ellsympow_goodred3(GEN E, GEN F, GEN p, long m, long vN, long *cnd, long *w)
1748 : {
1749 70 : long o = ellsympow_inertia3(E, vN);
1750 70 : long bet = ellsympow_betam(o, m);
1751 70 : *cnd = m + 1 - bet + ellsympow_deltam3(o, m, vN);
1752 70 : *w = odd(m)? ellsympow_rootno3(E, p, o, m): 1;
1753 70 : if (o==1 || o==2)
1754 0 : return ellsympow_abelian(p, ellap(F, p), m, o);
1755 70 : if ((o==3 || o==6) && ellsympow_isabelian3(F))
1756 0 : return ellsympow_abelian(p, p, m, o);
1757 : else
1758 70 : return ellsympow_nonabelian(p, m, bet);
1759 : }
1760 :
1761 : static long
1762 28 : ellsympow_inertia2(GEN F, long vN)
1763 : {
1764 28 : long vM = itos(gel(elllocalred(F, gen_2),1));
1765 28 : GEN c6 = ell_get_c6(F);
1766 28 : long v6 = signe(c6) ? vali(c6): 24;
1767 28 : if (vM==0) return vN==0 ? 1: 2;
1768 28 : if (vM==2) return vN==2 ? 3: 6;
1769 14 : if (vM==5) return 8;
1770 7 : if (vM==8) return v6>=9? 8: 4;
1771 7 : if (vM==3 || vN==7) return 24;
1772 0 : return 0;
1773 : }
1774 :
1775 : static long
1776 28 : ellsympow_deltam2(long o, long m, long vN)
1777 : {
1778 28 : if ((o==2 || o==6) && vN==4) return ellsympow_epsm(2, m);
1779 28 : if ((o==2 || o==6) && vN==6) return 2*ellsympow_epsm(2, m);
1780 28 : if (o==4) return 2*ellsympow_epsm(4, m)+ellsympow_epsm(2, m);
1781 28 : if (o==8 && vN==5) return ellsympow_epsm(8, m)+ellsympow_epsm(2, m)/2;
1782 21 : if (o==8 && vN==6) return ellsympow_epsm(8, m)+ellsympow_epsm(2, m);
1783 21 : if (o==8 && vN==8) return ellsympow_epsm(8, m)+ellsympow_epsm(4, m)+ellsympow_epsm(2, m);
1784 21 : if (o==24 && vN==3) return (2*ellsympow_epsm(8, m)+ellsympow_epsm(2, m))/6;
1785 14 : if (o==24 && vN==4) return (ellsympow_epsm(8, m)+ellsympow_epsm(2, m)*2)/3;
1786 14 : if (o==24 && vN==6) return (ellsympow_epsm(8, m)+ellsympow_epsm(2, m)*5)/3;
1787 14 : if (o==24 && vN==7) return (ellsympow_epsm(8, m)*10+ellsympow_epsm(2, m)*5)/6;
1788 14 : return 0;
1789 : }
1790 :
1791 : static long
1792 0 : ellsympow_isabelian2(GEN F)
1793 0 : { return umodi2n(ell_get_c4(F),7) == 96; }
1794 :
1795 : static long
1796 0 : ellsympow_rootno2(GEN E, long vN, long m, long bet)
1797 : {
1798 0 : long eps2 = (m + 1 - bet)>>1;
1799 0 : long eta = odd(vN) && m%8==3 ? -1 : 1;
1800 0 : long w2 = odd(eps2) ? ellrootno(E, gen_2): 1;
1801 0 : return eta == w2 ? 1 : -1;
1802 : }
1803 :
1804 : static GEN
1805 28 : ellsympow_goodred2(GEN E, GEN F, GEN p, long m, long vN, long *cnd, long *w)
1806 : {
1807 28 : long o = ellsympow_inertia2(F, vN);
1808 28 : long bet = ellsympow_betam(o, m);
1809 28 : *cnd = m + 1 - bet + ellsympow_deltam2(o, m, vN);
1810 28 : *w = odd(m) ? ellsympow_rootno2(E, vN, m, bet): 1;
1811 28 : if (o==1 || o==2)
1812 0 : return ellsympow_abelian(p, ellap(F, p), m, o);
1813 28 : if (o==4 && ellsympow_isabelian2(F))
1814 0 : return ellsympow_abelian(p, p, m, o);
1815 : else
1816 28 : return ellsympow_nonabelian(p, m, bet);
1817 : }
1818 :
1819 : static GEN
1820 189 : ellminimaldotwist(GEN E, GEN *pD)
1821 : {
1822 189 : GEN D = ellminimaltwistcond(E), Et = elltwist(E, D), Etmin;
1823 189 : if (pD) *pD = D;
1824 189 : Etmin = ellminimalmodel(Et, NULL);
1825 189 : obj_free(Et); return Etmin;
1826 : }
1827 :
1828 : /* Based on
1829 : Symmetric powers of elliptic curve L-functions,
1830 : Phil Martin and Mark Watkins, ANTS VII
1831 : <http://magma.maths.usyd.edu.au/users/watkins/papers/antsVII.pdf>
1832 : with thanks to Mark Watkins. BA20180402
1833 : */
1834 : static GEN
1835 140 : lfunellsympow(GEN e, ulong m)
1836 : {
1837 140 : pari_sp av = avma;
1838 : GEN B, N, Nfa, pr, ex, ld, bad, ejd, et, pole;
1839 140 : long i, l, mero, w = (m&7)==1 || (m&7)==3 ? -1: 1;
1840 140 : checkell_Q(e);
1841 140 : e = ellminimalmodel(e, NULL);
1842 140 : ejd = Q_denom(ell_get_j(e));
1843 140 : mero = m==0 || (m%4==0 && ellQ_get_CM(e)<0);
1844 140 : ellQ_get_Nfa(e, &N, &Nfa);
1845 140 : pr = gel(Nfa,1);
1846 140 : ex = gel(Nfa,2); l = lg(pr);
1847 140 : if (ugcd(umodiu(N,6), 6) == 1)
1848 21 : et = NULL;
1849 : else
1850 119 : et = ellminimaldotwist(e, NULL);
1851 140 : B = gen_1;
1852 140 : bad = cgetg(l, t_VEC);
1853 336 : for (i=1; i<l; i++)
1854 : {
1855 196 : long vN = itos(gel(ex,i));
1856 196 : GEN p = gel(pr,i), eul;
1857 : long cnd, wp;
1858 196 : if (dvdii(ejd, p))
1859 98 : eul = ellsympow_multred(e, p, m, vN, &cnd, &wp);
1860 98 : else if (equaliu(p, 2))
1861 28 : eul = ellsympow_goodred2(e, et, p, m, vN, &cnd, &wp);
1862 70 : else if (equaliu(p, 3))
1863 70 : eul = ellsympow_goodred3(e, et, p, m, vN, &cnd, &wp);
1864 : else
1865 0 : eul = ellsympow_goodred(e, p, m, &cnd, &wp);
1866 196 : gel(bad, i) = mkvec2(p, ginv(eul));
1867 196 : B = mulii(B, powiu(p,cnd));
1868 196 : w *= wp;
1869 : }
1870 140 : pole = mero ? mkvec(mkvec2(stoi(1+(m>>1)),gen_0)): NULL;
1871 280 : ld = mkvecn(mero? 7: 6, tag(mkvec2(mkvec2(e,utoi(m)),bad), t_LFUN_SYMPOW_ELL),
1872 140 : gen_0, ellsympow_gamma(m), stoi(m+1), B, stoi(w), pole);
1873 140 : if (et) obj_free(et);
1874 140 : return gerepilecopy(av, ld);
1875 : }
1876 :
1877 : GEN
1878 70 : lfunsympow(GEN ldata, ulong m)
1879 : {
1880 70 : ldata = lfunmisc_to_ldata_shallow(ldata);
1881 70 : if (ldata_get_type(ldata) != t_LFUN_ELL)
1882 0 : pari_err_IMPL("lfunsympow");
1883 70 : return lfunellsympow(gel(ldata_get_an(ldata), 2), m);
1884 : }
1885 :
1886 : static GEN
1887 28 : lfunmfspec_i(GEN lmisc, long bit)
1888 : {
1889 : GEN linit, ldataf, v, ve, vo, om, op, B, dom;
1890 : long k, k2, j;
1891 :
1892 28 : ldataf = lfunmisc_to_ldata_shallow(lmisc);
1893 28 : if (!gequal(ldata_get_gammavec(ldataf), mkvec2(gen_0,gen_1)))
1894 0 : pari_err_TYPE("lfunmfspec", lmisc);
1895 28 : k = gtos(ldata_get_k(ldataf));
1896 28 : if (k == 1) return mkvec2(cgetg(1, t_VEC), gen_1);
1897 21 : dom = mkvec3(dbltor(k/2.), dbltor((k-2)/2.), gen_0);
1898 21 : if (is_linit(lmisc) && linit_get_type(lmisc) == t_LDESC_INIT
1899 0 : && sdomain_isincl((double)k, dom, lfun_get_dom(linit_get_tech(lmisc))))
1900 0 : linit = lmisc;
1901 : else
1902 21 : linit = lfuninit(ldataf, dom, 0, bit);
1903 21 : B = int2n(bit/4);
1904 21 : v = cgetg(k, t_VEC);
1905 168 : for (j = 1; j < k; j++) gel(v,j) = lfunlambda(linit, utoi(j), bit);
1906 21 : om = gel(v,1);
1907 21 : if (odd(k)) return mkvec2(bestappr(gdiv(v, om), B), om);
1908 :
1909 7 : k2 = k/2;
1910 7 : ve = cgetg(k2, t_VEC);
1911 7 : vo = cgetg(k2+1, t_VEC);
1912 7 : gel(vo,1) = om;
1913 42 : for (j = 1; j < k2; j++)
1914 : {
1915 35 : gel(ve,j) = gel(v,2*j);
1916 35 : gel(vo,j+1) = gel(v,2*j+1);
1917 : }
1918 7 : if (k2 == 1) { om = gen_1; op = gel(v,1); }
1919 7 : else { om = gel(v,2); op = gel(v,3); }
1920 7 : if (maxss(gexpo(imag_i(om)), gexpo(imag_i(op))) > -bit/2)
1921 0 : pari_err_TYPE("lfunmfspec", lmisc);
1922 7 : ve = gdiv(ve, om);
1923 7 : vo = gdiv(vo, op);
1924 7 : return mkvec4(bestappr(ve,B), bestappr(vo,B), om, op);
1925 : }
1926 : GEN
1927 28 : lfunmfspec(GEN lmisc, long bit)
1928 : {
1929 28 : pari_sp av = avma;
1930 28 : return gerepilecopy(av, lfunmfspec_i(lmisc, bit));
1931 : }
1932 :
1933 : static long
1934 28 : ellsymsq_bad2(GEN c4, GEN c6, long e)
1935 : {
1936 28 : switch (e)
1937 : {
1938 14 : case 2: return 1;
1939 7 : case 3: return 0;
1940 7 : case 5: return 0;
1941 0 : case 7: return 0;
1942 0 : case 8:
1943 0 : if (!umodi2n(c6,9)) return 0;
1944 0 : return umodi2n(c4,7)==32 ? 1 : -1;
1945 0 : default: return 0;
1946 : }
1947 : }
1948 : static long
1949 14 : ellsymsq_bad3(GEN c4, GEN c6, long e)
1950 : {
1951 : long c6_243, c4_81;
1952 14 : switch (e)
1953 : {
1954 0 : case 2: return 1;
1955 14 : case 3: return 0;
1956 0 : case 5: return 0;
1957 0 : case 4:
1958 0 : c4_81 = umodiu(c4,81);
1959 0 : if (c4_81 == 27) return -1;
1960 0 : if (c4_81%27 != 9) return 1;
1961 0 : c6_243 = umodiu(c6,243);
1962 0 : return (c6_243==108 || c6_243==135)? -1: 1;
1963 0 : default: return 0;
1964 : }
1965 : }
1966 : static int
1967 0 : c4c6_testp(GEN c4, GEN c6, GEN p)
1968 0 : { GEN p2 = sqri(p); return (dvdii(c6,p2) && !dvdii(c4,p2)); }
1969 : /* assume e = v_p(N) >= 2 */
1970 : static long
1971 42 : ellsymsq_badp(GEN c4, GEN c6, GEN p, long e)
1972 : {
1973 42 : if (absequaliu(p, 2)) return ellsymsq_bad2(c4, c6, e);
1974 14 : if (absequaliu(p, 3)) return ellsymsq_bad3(c4, c6, e);
1975 0 : switch(umodiu(p, 12UL))
1976 : {
1977 0 : case 1: return -1;
1978 0 : case 5: return c4c6_testp(c4,c6,p)? -1: 1;
1979 0 : case 7: return c4c6_testp(c4,c6,p)? 1:-1;
1980 0 : default:return 1; /* p%12 = 11 */
1981 : }
1982 : }
1983 : static GEN
1984 70 : lfunellsymsqmintwist(GEN e)
1985 : {
1986 70 : pari_sp av = avma;
1987 : GEN N, Nfa, P, E, V, c4, c6, ld;
1988 : long i, l, k;
1989 70 : checkell_Q(e);
1990 70 : e = ellminimalmodel(e, NULL);
1991 70 : ellQ_get_Nfa(e, &N, &Nfa);
1992 70 : c4 = ell_get_c4(e);
1993 70 : c6 = ell_get_c6(e);
1994 70 : P = gel(Nfa,1); l = lg(P);
1995 70 : E = gel(Nfa,2);
1996 70 : V = cgetg(l, t_VEC);
1997 196 : for (i=k=1; i<l; i++)
1998 : {
1999 126 : GEN p = gel(P,i);
2000 126 : long a, e = itos(gel(E,i));
2001 126 : if (e == 1) continue;
2002 42 : a = ellsymsq_badp(c4, c6, p, e);
2003 42 : gel(V,k++) = mkvec2(p, stoi(a));
2004 : }
2005 70 : setlg(V, k);
2006 70 : ld = lfunellsympow(e, 2);
2007 70 : return gerepilecopy(av, mkvec2(ld, V));
2008 : }
2009 :
2010 : static GEN
2011 70 : mfpeters(GEN ldata2, GEN fudge, GEN N, long k, long bitprec)
2012 : {
2013 70 : GEN t, L = real_i(lfun(ldata2, stoi(k), bitprec));
2014 70 : long prec = nbits2prec(bitprec);
2015 70 : t = powrs(mppi(prec), k+1); shiftr_inplace(t, 2*k-1); /* Pi/2 * (4Pi)^k */
2016 70 : return gmul(gdiv(gmul(mulii(N,mpfact(k-1)), fudge), t), L);
2017 : }
2018 :
2019 : /* Assume E to be twist-minimal */
2020 : static GEN
2021 70 : lfunellmfpetersmintwist(GEN E, long bitprec)
2022 : {
2023 70 : pari_sp av = avma;
2024 70 : GEN symsq, veceuler, N = ellQ_get_N(E), fudge = gen_1;
2025 70 : long j, k = 2;
2026 70 : symsq = lfunellsymsqmintwist(E);
2027 70 : veceuler = gel(symsq,2);
2028 112 : for (j = 1; j < lg(veceuler); j++)
2029 : {
2030 42 : GEN v = gel(veceuler,j), p = gel(v,1), q = powis(p,1-k);
2031 42 : long s = signe(gel(v,2));
2032 42 : if (s) fudge = gmul(fudge, s==1 ? gaddsg(1, q): gsubsg(1, q));
2033 : }
2034 70 : return gerepileupto(av, mfpeters(gel(symsq,1),fudge,N,k,bitprec));
2035 : }
2036 :
2037 : /* From Christophe Delaunay, http://delaunay.perso.math.cnrs.fr/these.pdf */
2038 : static GEN
2039 70 : elldiscfix(GEN E, GEN Et, GEN D)
2040 : {
2041 70 : GEN N = ellQ_get_N(E), Nt = ellQ_get_N(Et);
2042 70 : GEN P = gel(absZ_factor(D), 1);
2043 70 : GEN f = gen_1;
2044 70 : long i, l = lg(P);
2045 133 : for (i=1; i < l; i++)
2046 : {
2047 63 : GEN r, p = gel(P,i);
2048 63 : long v = Z_pval(N, p), vt = Z_pval(Nt, p);
2049 63 : if (v <= vt) continue;
2050 : /* v > vt */
2051 49 : if (absequaliu(p, 2))
2052 : {
2053 28 : if (vt == 0 && v >= 4)
2054 0 : r = shifti(subsi(9, sqri(ellap(Et, p))), v-3); /* 9=(2+1)^2 */
2055 28 : else if (vt == 1)
2056 7 : r = gmul2n(utoipos(3), v-3); /* not in Z if v=2 */
2057 21 : else if (vt >= 2)
2058 21 : r = int2n(v-vt);
2059 : else
2060 0 : r = gen_1; /* vt = 0, 1 <= v <= 3 */
2061 : }
2062 21 : else if (vt >= 1)
2063 14 : r = gdiv(subiu(sqri(p), 1), p);
2064 : else
2065 7 : r = gdiv(mulii(subiu(p, 1), subii(sqri(addiu(p, 1)), sqri(ellap(Et, p)))), p);
2066 49 : f = gmul(f, r);
2067 : }
2068 70 : return f;
2069 : }
2070 :
2071 : GEN
2072 70 : lfunellmfpeters(GEN E, long bitprec)
2073 : {
2074 70 : pari_sp ltop = avma;
2075 70 : GEN D, Et = ellminimaldotwist(E, &D);
2076 70 : GEN nor = lfunellmfpetersmintwist(Et, bitprec);
2077 70 : GEN nor2 = gmul(nor, elldiscfix(E, Et, D));
2078 70 : obj_free(Et); return gerepileupto(ltop, nor2);
2079 : }
2080 :
2081 : /*************************************************************/
2082 : /* Genus 2 curves */
2083 : /*************************************************************/
2084 :
2085 : static void
2086 233611 : Flv_diffnext(GEN d, ulong p)
2087 : {
2088 233611 : long j, n = lg(d)-1;
2089 1635277 : for(j = n; j>=2; j--)
2090 1401666 : uel(d,j) = Fl_add(uel(d,j), uel(d,j-1), p);
2091 233611 : }
2092 :
2093 : static GEN
2094 1995 : Flx_difftable(GEN P, ulong p)
2095 : {
2096 1995 : long i, n = degpol(P);
2097 1995 : GEN V = cgetg(n+2, t_VECSMALL);
2098 1995 : uel(V, n+1) = Flx_constant(P);
2099 13965 : for(i = n; i >= 1; i--)
2100 : {
2101 11970 : P = Flx_diff1(P, p);
2102 11970 : uel(V, i) = Flx_constant(P);
2103 : }
2104 1995 : return V;
2105 : }
2106 :
2107 : static long
2108 1995 : Flx_genus2trace_naive(GEN H, ulong p)
2109 : {
2110 1995 : pari_sp av = avma;
2111 : ulong i, j;
2112 1995 : long a, n = degpol(H);
2113 1995 : GEN k = const_vecsmall(p, -1), d;
2114 1995 : k[1] = 0;
2115 117803 : for (i=1, j=1; i < p; i += 2, j = Fl_add(j, i, p))
2116 115808 : k[j+1] = 1;
2117 1995 : a = n == 5 ? 0: k[1+Flx_lead(H)];
2118 1995 : d = Flx_difftable(H, p);
2119 235606 : for (i=0; i < p; i++)
2120 : {
2121 233611 : a += k[1+uel(d,n+1)];
2122 233611 : Flv_diffnext(d, p);
2123 : }
2124 1995 : set_avma(av);
2125 1995 : return a;
2126 : }
2127 :
2128 : static GEN
2129 2156 : dirgenus2(GEN Q, GEN p, long n)
2130 : {
2131 2156 : pari_sp av = avma;
2132 : GEN f;
2133 2156 : if (n > 2)
2134 161 : f = RgX_recip(hyperellcharpoly(gmul(Q,gmodulo(gen_1, p))));
2135 : else
2136 : {
2137 1995 : ulong pp = itou(p);
2138 1995 : GEN Qp = ZX_to_Flx(Q, pp);
2139 1995 : long t = Flx_genus2trace_naive(Qp, pp);
2140 1995 : f = deg1pol_shallow(stoi(t), gen_1, 0);
2141 : }
2142 2156 : return gerepileupto(av, RgXn_inv_i(f, n));
2143 : }
2144 :
2145 : GEN
2146 875 : dirgenus2_worker(GEN P, ulong X, GEN Q)
2147 : {
2148 875 : pari_sp av = avma;
2149 875 : long i, l = lg(P);
2150 875 : GEN V = cgetg(l, t_VEC);
2151 3031 : for(i = 1; i < l; i++)
2152 : {
2153 2156 : ulong p = uel(P,i);
2154 2156 : long d = ulogint(X, p) + 1; /* minimal d such that p^d > X */
2155 2156 : gel(V,i) = dirgenus2(Q, utoi(uel(P,i)), d);
2156 : }
2157 875 : return gerepilecopy(av, mkvec2(P,V));
2158 : }
2159 :
2160 : static GEN
2161 196 : vecan_genus2(GEN an, long L)
2162 : {
2163 196 : GEN Q = gel(an,1), bad = gel(an, 2);
2164 196 : GEN worker = snm_closure(is_entry("_dirgenus2_worker"), mkvec(Q));
2165 196 : return pardireuler(worker, gen_2, stoi(L), NULL, bad);
2166 : }
2167 :
2168 : static GEN
2169 0 : eulerf_genus2(GEN an, GEN p)
2170 : {
2171 0 : GEN Q = gel(an,1), bad = gel(an, 2);
2172 0 : GEN f = eulerf_bad(bad, p);
2173 0 : if (f) return f;
2174 0 : f = RgX_recip(hyperellcharpoly(gmul(Q,gmodulo(gen_1, p))));
2175 0 : return mkrfrac(gen_1,f);
2176 : }
2177 :
2178 : static GEN
2179 49 : genus2_redmodel(GEN P, GEN p)
2180 : {
2181 49 : GEN M = FpX_factor(P, p);
2182 49 : GEN F = gel(M,1), E = gel(M,2);
2183 49 : long i, k, r = lg(F);
2184 49 : GEN U = scalarpol(leading_coeff(P), varn(P));
2185 49 : GEN G = cgetg(r, t_COL);
2186 161 : for (i=1, k=0; i<r; i++)
2187 : {
2188 112 : if (E[i]>1)
2189 56 : gel(G,++k) = gel(F,i);
2190 112 : if (odd(E[i]))
2191 77 : U = FpX_mul(U, gel(F,i), p);
2192 : }
2193 49 : setlg(G,++k);
2194 49 : return mkvec2(G,U);
2195 : }
2196 :
2197 : static GEN
2198 308 : oneminusxd(long d)
2199 : {
2200 308 : return gsub(gen_1, pol_xn(d, 0));
2201 : }
2202 :
2203 : static GEN
2204 21 : ellfromeqncharpoly(GEN P, GEN Q, GEN p)
2205 : {
2206 : long v;
2207 : GEN E, F, t, y;
2208 21 : v = fetch_var();
2209 21 : y = pol_x(v);
2210 21 : F = gsub(gadd(ZX_sqr(y), gmul(y, Q)), P);
2211 21 : E = ellinit(ellfromeqn(F), p, DEFAULTPREC);
2212 21 : delete_var();
2213 21 : t = ellap(E, p);
2214 21 : obj_free(E);
2215 21 : return mkpoln(3, gen_1, negi(t), p);
2216 : }
2217 :
2218 : static GEN
2219 0 : nfellcharpoly(GEN e, GEN T, GEN p)
2220 : {
2221 : GEN nf, E, t;
2222 0 : nf = nfinit(mkvec2(T, mkvec(p)), DEFAULTPREC);
2223 0 : E = ellinit(e, nf, DEFAULTPREC);
2224 0 : if (lg(E)==1) return NULL;
2225 0 : t = elleulerf(E, p);
2226 0 : obj_free(E);
2227 0 : return t;
2228 : }
2229 :
2230 : static GEN
2231 0 : genus2_red5(GEN P, GEN T, GEN p)
2232 : {
2233 0 : long vx = varn(P), vy = varn(T);
2234 0 : GEN f = shallowcopy(T);
2235 0 : setvarn(f, vx);
2236 : while(1)
2237 0 : {
2238 : GEN Pr, R, r, Rs;
2239 0 : (void) ZXX_pvalrem(P, p, &Pr);
2240 0 : R = FpXQX_roots_mult(Pr, 3, T, p);
2241 0 : if (lg(R)==1) return P;
2242 0 : r = FpX_center(gel(R,1), p, shifti(p,-1));
2243 0 : r = gel(R,1);
2244 0 : P = RgX_affine(P, p, r);
2245 0 : setvarn(r, vx);
2246 0 : f = RgX_Rg_div(gsub(f, r), p);
2247 0 : Rs = RgX_rem(RgXY_swap(P, 3, vy), gsub(f, pol_x(vy)));
2248 0 : P = RgXY_swap(Rs, 3, vy);
2249 0 : (void) ZXX_pvalrem(P, sqri(p), &P);
2250 : }
2251 : }
2252 :
2253 : static GEN
2254 63 : genus2_type5(GEN P, GEN p)
2255 : {
2256 : GEN E, F, T, a, a2, Q;
2257 : long v;
2258 63 : (void) ZX_pvalrem(P, p, &F);
2259 63 : F = FpX_red(F, p);
2260 63 : if (degpol(F) < 1) return NULL;
2261 63 : F = FpX_factor(F, p);
2262 63 : if (mael(F,2,1) != 3 || degpol(gmael(F,1,1)) != 2) return NULL;
2263 0 : T = gmael(F, 1, 1);
2264 0 : v = fetch_var_higher();
2265 0 : Q = RgV_to_RgX(ZX_digits(P, T), v);
2266 0 : Q = genus2_red5(Q, T, p);
2267 0 : a = gel(Q,5); a2 = ZX_sqr(a);
2268 0 : E = mkvec5(gen_0, gel(Q,4), gen_0, ZX_mul(gel(Q,3),a), ZX_mul(gel(Q,2),a2));
2269 0 : delete_var();
2270 0 : return nfellcharpoly(E, T, p);
2271 : }
2272 :
2273 : /* Assume P has semistable reduction at p */
2274 : static GEN
2275 49 : genus2_eulerfact_semistable(GEN P, GEN p)
2276 : {
2277 49 : GEN Pp = FpX_red(P, p);
2278 49 : GEN GU = genus2_redmodel(Pp, p);
2279 49 : long d = 6-degpol(Pp), v = d/2, w = odd(d);
2280 : GEN abe, tor;
2281 49 : GEN ki, kp = pol_1(0), kq = pol_1(0);
2282 49 : GEN F = gel(GU,1), Q = gel(GU,2);
2283 49 : long dQ = degpol(Q), lF = lg(F)-1;
2284 :
2285 0 : abe = dQ >= 5 ? RgX_recip(hyperellcharpoly(gmul(Q,gmodulo(gen_1,p))))
2286 98 : : dQ >= 3 ? RgX_recip(ellfromeqncharpoly(Q,gen_0,p))
2287 49 : : pol_1(0);
2288 42 : ki = dQ != 0 ? oneminusxd(1)
2289 56 : : Fp_issquare(gel(Q,2),p) ? ZX_sqr(oneminusxd(1))
2290 7 : : oneminusxd(2);
2291 49 : if (lF)
2292 : {
2293 : long i;
2294 98 : for(i=1; i <= lF; i++)
2295 : {
2296 56 : GEN Fi = gel(F, i);
2297 56 : long d = degpol(Fi);
2298 56 : GEN e = FpX_rem(Q, Fi, p);
2299 91 : GEN kqf = lgpol(e)==0 ? oneminusxd(d):
2300 70 : FpXQ_issquare(e, Fi, p) ? ZX_sqr(oneminusxd(d))
2301 35 : : oneminusxd(2*d);
2302 56 : kp = gmul(kp, oneminusxd(d));
2303 56 : kq = gmul(kq, kqf);
2304 : }
2305 : }
2306 49 : if (v)
2307 : {
2308 7 : GEN kqoo = w==1 ? oneminusxd(1):
2309 0 : Fp_issquare(leading_coeff(Q), p)? ZX_sqr(oneminusxd(1))
2310 0 : : oneminusxd(2);
2311 7 : kp = gmul(kp, oneminusxd(1));
2312 7 : kq = gmul(kq, kqoo);
2313 : }
2314 49 : tor = RgX_div(ZX_mul(oneminusxd(1), kq), ZX_mul(ki, kp));
2315 49 : return ginv( ZX_mul(abe, tor) );
2316 : }
2317 :
2318 : static GEN
2319 49 : genus2_eulerfact(GEN P, GEN p)
2320 : {
2321 49 : pari_sp av = avma;
2322 49 : GEN W, R = genus2_type5(P, p), E;
2323 49 : if (R) return R;
2324 49 : W = hyperellextremalmodels(P, 2, p);
2325 49 : if (lg(W) < 3) return genus2_eulerfact_semistable(P,p);
2326 0 : E = gmul(genus2_eulerfact_semistable(gel(W,1),p),
2327 0 : genus2_eulerfact_semistable(gel(W,2),p));
2328 0 : return gerepileupto(av, E);
2329 : }
2330 :
2331 : static GEN
2332 28 : F2x_genus2_find_trans(GEN P, GEN Q, GEN F)
2333 : {
2334 28 : pari_sp av = avma;
2335 28 : long i, d = F2x_degree(F), v = P[1];
2336 : GEN M, C, V;
2337 28 : M = cgetg(d+1, t_MAT);
2338 84 : for (i=1; i<=d; i++)
2339 : {
2340 56 : GEN Mi = F2x_rem(F2x_add(F2x_shift(Q,i-1), monomial_F2x(2*i-2,v)), F);
2341 56 : gel(M,i) = F2x_to_F2v(Mi, d);
2342 : }
2343 28 : C = F2x_to_F2v(F2x_rem(P, F), d);
2344 28 : V = F2m_F2c_invimage(M, C);
2345 28 : return gerepileuptoleaf(av, F2v_to_F2x(V, v));
2346 : }
2347 :
2348 : static GEN
2349 35 : F2x_genus2_trans(GEN P, GEN Q, GEN H)
2350 : {
2351 35 : return F2x_add(P,F2x_add(F2x_mul(H,Q), F2x_sqr(H)));
2352 : }
2353 :
2354 : static GEN
2355 42 : F2x_genus_redoo(GEN P, GEN Q, long k)
2356 : {
2357 42 : if (F2x_degree(P)==2*k)
2358 : {
2359 14 : long c = F2x_coeff(P,2*k-1), dQ = F2x_degree(Q);
2360 14 : if ((dQ==k-1 && c==1) || (dQ<k-1 && c==0))
2361 7 : return F2x_genus2_trans(P, Q, monomial_F2x(k, P[1]));
2362 : }
2363 35 : return P;
2364 : }
2365 :
2366 : static GEN
2367 35 : F2x_pseudodisc(GEN P, GEN Q)
2368 : {
2369 35 : GEN dP = F2x_deriv(P), dQ = F2x_deriv(Q);
2370 35 : return F2x_gcd(Q, F2x_add(F2x_mul(P, F2x_sqr(dQ)), F2x_sqr(dP)));
2371 : }
2372 :
2373 : static GEN
2374 14 : F2x_genus_red(GEN P, GEN Q)
2375 : {
2376 : long dP, dQ;
2377 : GEN F, FF;
2378 14 : P = F2x_genus_redoo(P, Q, 3);
2379 14 : P = F2x_genus_redoo(P, Q, 2);
2380 14 : P = F2x_genus_redoo(P, Q, 1);
2381 14 : dP = F2x_degree(P);
2382 14 : dQ = F2x_degree(Q);
2383 14 : FF = F = F2x_pseudodisc(P,Q);
2384 35 : while(F2x_degree(F)>0)
2385 : {
2386 21 : GEN M = gel(F2x_factor(F),1);
2387 21 : long i, l = lg(M);
2388 49 : for(i=1; i<l; i++)
2389 : {
2390 28 : GEN R = F2x_sqr(gel(M,i));
2391 28 : GEN H = F2x_genus2_find_trans(P, Q, R);
2392 28 : P = F2x_div(F2x_genus2_trans(P, Q, H), R);
2393 28 : Q = F2x_div(Q, gel(M,i));
2394 : }
2395 21 : F = F2x_pseudodisc(P, Q);
2396 : }
2397 14 : return mkvec4(P,Q,FF,mkvecsmall2(dP,dQ));
2398 : }
2399 :
2400 : /* Number of solutions of x^2+b*x+c */
2401 : static long
2402 21 : F2xqX_quad_nbroots(GEN b, GEN c, GEN T)
2403 : {
2404 21 : if (lgpol(b) > 0)
2405 : {
2406 14 : GEN d = F2xq_div(c, F2xq_sqr(b, T), T);
2407 14 : return F2xq_trace(d, T)? 0: 2;
2408 : }
2409 : else
2410 7 : return 1;
2411 : }
2412 :
2413 : static GEN
2414 14 : genus2_eulerfact2_semistable(GEN PQ)
2415 : {
2416 14 : GEN V = F2x_genus_red(ZX_to_F2x(gel(PQ, 1)), ZX_to_F2x(gel(PQ, 2)));
2417 14 : GEN P = gel(V, 1), Q = gel(V, 2);
2418 14 : GEN F = gel(V, 3), v = gel(V, 4);
2419 : GEN abe, tor;
2420 14 : GEN ki, kp = pol_1(0), kq = pol_1(0);
2421 14 : long dP = F2x_degree(P), dQ = F2x_degree(Q), d = maxss(dP, 2*dQ);
2422 14 : if (!lgpol(F)) return pol_1(0);
2423 21 : ki = dQ!=0 || dP>0 ? oneminusxd(1):
2424 7 : dP==-1 ? ZX_sqr(oneminusxd(1)): oneminusxd(2);
2425 28 : abe = d>=5? RgX_recip(hyperellcharpoly(gmul(PQ,gmodulss(1,2)))):
2426 14 : d>=3? RgX_recip(ellfromeqncharpoly(F2x_to_ZX(P), F2x_to_ZX(Q), gen_2)):
2427 14 : pol_1(0);
2428 14 : if (lgpol(F))
2429 : {
2430 14 : GEN M = gel(F2x_factor(F), 1);
2431 14 : long i, lF = lg(M)-1;
2432 35 : for(i=1; i <= lF; i++)
2433 : {
2434 21 : GEN Fi = gel(M, i);
2435 21 : long d = F2x_degree(Fi);
2436 21 : long nb = F2xqX_quad_nbroots(F2x_rem(Q, Fi), F2x_rem(P, Fi), Fi);
2437 35 : GEN kqf = nb==1 ? oneminusxd(d):
2438 0 : nb==2 ? ZX_sqr(oneminusxd(d))
2439 14 : : oneminusxd(2*d);
2440 21 : kp = gmul(kp, oneminusxd(d));
2441 21 : kq = gmul(kq, kqf);
2442 : }
2443 : }
2444 14 : if (maxss(v[1],2*v[2])<5)
2445 : {
2446 14 : GEN kqoo = v[1]>2*v[2] ? oneminusxd(1):
2447 0 : v[1]<2*v[2] ? ZX_sqr(oneminusxd(1))
2448 7 : : oneminusxd(2);
2449 7 : kp = gmul(kp, oneminusxd(1));
2450 7 : kq = gmul(kq, kqoo);
2451 : }
2452 14 : tor = RgX_div(ZX_mul(oneminusxd(1),kq), ZX_mul(ki, kp));
2453 14 : return ginv( ZX_mul(abe, tor) );
2454 : }
2455 :
2456 : static GEN
2457 14 : genus2_eulerfact2(GEN F, GEN PQ)
2458 : {
2459 14 : pari_sp av = avma;
2460 14 : GEN W, R = genus2_type5(F, gen_2), E;
2461 14 : if (R) return R;
2462 14 : W = hyperellextremalmodels(PQ, 2, gen_2);
2463 14 : if (lg(W) < 3) return genus2_eulerfact2_semistable(PQ);
2464 0 : E = gmul(genus2_eulerfact2_semistable(gel(W,1)),
2465 0 : genus2_eulerfact2_semistable(gel(W,2)));
2466 0 : return gerepileupto(av, E);
2467 : }
2468 :
2469 : GEN
2470 35 : lfungenus2(GEN G)
2471 : {
2472 35 : pari_sp ltop = avma;
2473 : GEN Ldata;
2474 35 : GEN gr = genus2red(G, NULL);
2475 35 : GEN N = gel(gr, 1), M = gel(gr, 2), PQ = gel(gr, 3), L = gel(gr, 4);
2476 35 : GEN e, F = gadd(gsqr(gel(PQ, 2)), gmul2n(gel(PQ, 1), 2));
2477 35 : long i, lL = lg(L), ram2;
2478 35 : ram2 = absequaliu(gmael(M,1,1),2);
2479 35 : if (ram2 && equalis(gmael(M,2,1),-1))
2480 14 : pari_warn(warner,"unknown valuation of conductor at 2");
2481 35 : e = cgetg(lL+(ram2?0:1), t_VEC);
2482 56 : gel(e,1) = mkvec2(gen_2, ram2 ? genus2_eulerfact2(F, PQ)
2483 21 : : ginv( RgX_recip(hyperellcharpoly(gmul(PQ,gmodulss(1,2))))) );
2484 84 : for(i = ram2? 2: 1; i < lL; i++)
2485 : {
2486 49 : GEN Li = gel(L, i);
2487 49 : GEN p = gel(Li, 1);
2488 49 : gel(e, ram2 ? i: i+1) = mkvec2(p, genus2_eulerfact(F,p));
2489 : }
2490 35 : Ldata = mkvecn(6, tag(mkvec2(F,e), t_LFUN_GENUS2),
2491 : gen_0, mkvec4(gen_0, gen_0, gen_1, gen_1), gen_2, N, gen_0);
2492 35 : return gerepilecopy(ltop, Ldata);
2493 : }
2494 :
2495 : /*************************************************************/
2496 : /* ETA QUOTIENTS */
2497 : /* An eta quotient is a matrix with 2 columns [m, r_m] with */
2498 : /* m >= 1 representing f(\tau)=\prod_m\eta(m\tau)^{r_m}. */
2499 : /*************************************************************/
2500 :
2501 : /* eta(x^v) + O(x^L) */
2502 : GEN
2503 1222 : eta_ZXn(long v, long L)
2504 : {
2505 1222 : long n, k = 0, v2 = 2*v, bn = v, cn = 0;
2506 : GEN P;
2507 1222 : if (!L) return zeropol(0);
2508 1222 : P = cgetg(L+2,t_POL); P[1] = evalsigne(1);
2509 72543 : for(n = 0; n < L; n++) gel(P,n+2) = gen_0;
2510 1222 : for(n = 0;; n++, bn += v2, cn += v)
2511 3000 : { /* k = v * (3*n-1) * n / 2; bn = v * (2*n+1); cn = v * n */
2512 : long k2;
2513 4222 : gel(P, k+2) = odd(n)? gen_m1: gen_1;
2514 4222 : k2 = k+cn; if (k2 >= L) break;
2515 3770 : k = k2;
2516 : /* k = v * (3*n+1) * n / 2 */;
2517 3770 : gel(P, k+2) = odd(n)? gen_m1: gen_1;
2518 3770 : k2 = k+bn; if (k2 >= L) break;
2519 3000 : k = k2;
2520 : }
2521 1222 : setlg(P, k+3); return P;
2522 : }
2523 : GEN
2524 322 : eta_product_ZXn(GEN eta, long L)
2525 : {
2526 322 : pari_sp av = avma;
2527 322 : GEN P = NULL, D = gel(eta,1), R = gel(eta,2);
2528 322 : long i, l = lg(D);
2529 1148 : for (i = 1; i < l; ++i)
2530 : {
2531 826 : GEN Q = eta_ZXn(D[i], L);
2532 826 : long r = R[i];
2533 826 : if (r < 0) { Q = RgXn_inv_i(Q, L); r = -r; }
2534 826 : if (r != 1) Q = RgXn_powu_i(Q, r, L);
2535 826 : P = P? ZXn_mul(P, Q, L): Q;
2536 826 : if (gc_needed(av,1) && i > 1)
2537 : {
2538 0 : if (DEBUGMEM>1) pari_warn(warnmem,"eta_product_ZXn");
2539 0 : P = gerepilecopy(av, P);
2540 : }
2541 : }
2542 322 : return P;
2543 : }
2544 : static GEN
2545 147 : vecan_eta(GEN an, long L)
2546 : {
2547 147 : long v = itos(gel(an, 3));
2548 : GEN t;
2549 147 : if (v > L) return zerovec(L);
2550 140 : t = eta_product_ZXn(an, L - v);
2551 140 : if (v) t = RgX_shift_shallow(t, v);
2552 140 : return RgX_to_RgV(t, L);
2553 : }
2554 : /* return 1 if cuspidal, 0 if holomorphic, -1 otherwise */
2555 : static int
2556 231 : etacuspidal(GEN N, GEN k, GEN B, GEN R, GEN NB)
2557 : {
2558 231 : long i, j, lD, l, cusp = 1;
2559 231 : pari_sp av = avma;
2560 : GEN D;
2561 231 : if (gsigne(k) < 0) return -1;
2562 224 : D = divisors(N); lD = lg(D); l = lg(B);
2563 1505 : for (i = 1; i < lD; i++)
2564 : {
2565 1281 : GEN t = gen_0, d = gel(D,i);
2566 : long s;
2567 3997 : for (j = 1; j < l; j++)
2568 2716 : t = addii(t, mulii(gel(NB,j), mulii(gel(R,j), sqri(gcdii(d, gel(B,j))))));
2569 1281 : s = signe(t);
2570 1281 : if (s < 0) return -1;
2571 1281 : if (!s) cusp = 0;
2572 : }
2573 224 : return gc_bool(av, cusp);
2574 : }
2575 : /* u | 24, level N = u*N0, N0 = lcm(B), NB[i] = N0/B[i] */
2576 : static int
2577 49 : etaselfdual(GEN B, GEN R, GEN NB, ulong u)
2578 : {
2579 49 : pari_sp av = avma;
2580 49 : long i, l = lg(B);
2581 161 : for (i = 1; i < l; i++)
2582 : {
2583 112 : long j = ZV_search(B, muliu(gel(NB,i), u)); /* search for N / B[i] */
2584 112 : set_avma(av); if (!j || !equalii(gel(R,i),gel(R,j))) return 0;
2585 : }
2586 49 : return 1;
2587 : }
2588 : /* return Nebentypus of eta quotient, k2 = 2*k integral */
2589 : static GEN
2590 203 : etachar(GEN B, GEN R, GEN k2)
2591 : {
2592 203 : long i, l = lg(B);
2593 203 : GEN P = gen_1;
2594 546 : for (i = 1; i < l; ++i) if (mpodd(gel(R,i))) P = mulii(P, gel(B,i));
2595 203 : switch(Mod4(k2))
2596 : {
2597 133 : case 0: break;
2598 42 : case 2: P = negi(P); break;
2599 28 : default: P = shifti(P, 1); break;
2600 : }
2601 203 : return coredisc(P);
2602 : }
2603 : /* Return 0 if not on gamma_0(N). Sets conductor, modular weight, character,
2604 : * canonical matrix, v_q(eta), sd = 1 iff self-dual, cusp = 1 iff cuspidal
2605 : * [0 if holomorphic at all cusps, else -1] */
2606 : long
2607 259 : etaquotype(GEN *peta, GEN *pN, GEN *pk, GEN *CHI, long *pv, long *sd,
2608 : long *cusp)
2609 : {
2610 259 : GEN B, R, S, T, N, NB, eta = *peta;
2611 : long l, i, u, S24;
2612 :
2613 259 : if (lg(eta) != 3) pari_err_TYPE("lfunetaquo", eta);
2614 259 : switch(typ(eta))
2615 : {
2616 77 : case t_VEC: eta = mkmat2(mkcol(gel(eta,1)), mkcol(gel(eta,2))); break;
2617 182 : case t_MAT: break;
2618 0 : default: pari_err_TYPE("lfunetaquo", eta);
2619 : }
2620 259 : if (!RgV_is_ZVpos(gel(eta,1)) || !RgV_is_ZV(gel(eta,2)))
2621 0 : pari_err_TYPE("lfunetaquo", eta);
2622 259 : *peta = eta = famat_reduce(eta);
2623 259 : B = gel(eta,1); l = lg(B); /* sorted in increasing order */
2624 259 : R = gel(eta,2);
2625 259 : N = ZV_lcm(B); NB = cgetg(l, t_VEC);
2626 721 : for (i = 1; i < l; i++) gel(NB,i) = diviiexact(N, gel(B,i));
2627 259 : S = gen_0; T = gen_0; u = 0;
2628 721 : for (i = 1; i < l; ++i)
2629 : {
2630 462 : GEN b = gel(B,i), r = gel(R,i);
2631 462 : S = addii(S, mulii(b, r));
2632 462 : T = addii(T, r);
2633 462 : u += umodiu(r,24) * umodiu(gel(NB,i), 24);
2634 : }
2635 259 : S = divis_rem(S, 24, &S24);
2636 259 : if (S24) return 0; /* nonintegral valuation at oo */
2637 252 : u = 24 / ugcd(24, u % 24);
2638 252 : *pN = muliu(N, u); /* level */
2639 252 : *pk = gmul2n(T,-1); /* weight */
2640 252 : *pv = itos(S); /* valuation */
2641 252 : if (cusp) *cusp = etacuspidal(*pN, *pk, B, R, NB);
2642 252 : if (sd) *sd = etaselfdual(B, R, NB, u);
2643 252 : if (CHI) *CHI = etachar(B, R, T);
2644 252 : return 1;
2645 : }
2646 :
2647 : GEN
2648 49 : lfunetaquo(GEN eta0)
2649 : {
2650 49 : pari_sp ltop = avma;
2651 49 : GEN Ldata, N, BR, k, eta = eta0;
2652 : long v, sd, cusp;
2653 49 : if (!etaquotype(&eta, &N, &k, NULL, &v, &sd, &cusp))
2654 0 : pari_err_TYPE("lfunetaquo", eta0);
2655 49 : if (!cusp) pari_err_IMPL("noncuspidal eta quotient");
2656 49 : if (!sd) pari_err_IMPL("non self-dual eta quotient");
2657 49 : if (typ(k) != t_INT) pari_err_TYPE("lfunetaquo [nonintegral weight]", eta0);
2658 49 : BR = mkvec3(ZV_to_zv(gel(eta,1)), ZV_to_zv(gel(eta,2)), stoi(v - 1));
2659 49 : Ldata = mkvecn(6, tag(BR,t_LFUN_ETA), gen_0, mkvec2(gen_0,gen_1), k,N, gen_1);
2660 49 : return gerepilecopy(ltop, Ldata);
2661 : }
2662 :
2663 : static GEN
2664 399 : vecan_qf(GEN Q, long L)
2665 : {
2666 399 : GEN v, w = qfrep0(Q, utoi(L), 1);
2667 : long i;
2668 399 : v = cgetg(L+1, t_VEC);
2669 26698 : for (i = 1; i <= L; i++) gel(v,i) = utoi(2 * w[i]);
2670 399 : return v;
2671 : }
2672 :
2673 : long
2674 336 : qfiseven(GEN M)
2675 : {
2676 336 : long i, l = lg(M);
2677 784 : for (i=1; i<l; i++)
2678 679 : if (mpodd(gcoeff(M,i,i))) return 0;
2679 105 : return 1;
2680 : }
2681 :
2682 : GEN
2683 91 : lfunqf(GEN M, long prec)
2684 : {
2685 91 : pari_sp ltop = avma;
2686 : long n;
2687 : GEN k, D, d, Mi, Ldata, poles, eno, dual;
2688 :
2689 91 : if (typ(M) != t_MAT) pari_err_TYPE("lfunqf", M);
2690 91 : if (!RgM_is_ZM(M)) pari_err_TYPE("lfunqf [not integral]", M);
2691 91 : n = lg(M)-1;
2692 91 : k = uutoQ(n,2);
2693 91 : M = Q_primpart(M);
2694 91 : Mi = ZM_inv(M, &d); /* d M^(-1) */
2695 91 : if (!qfiseven(M)) { M = gmul2n(M, 1); d = shifti(d,1); }
2696 91 : if (!qfiseven(Mi)){ Mi= gmul2n(Mi,1); d = shifti(d,1); }
2697 : /* det(Mi) = d^n/det(M), D^2 = det(Mi)/det(M) */
2698 91 : D = gdiv(gpow(d,k,prec), ZM_det(M));
2699 91 : if (!issquareall(D, &eno)) eno = gsqrt(D, prec);
2700 91 : dual = gequal1(D) ? gen_0: tag(Mi, t_LFUN_QF);
2701 91 : poles = mkcol2(mkvec2(k, simple_pole(gmul2n(eno,1))),
2702 : mkvec2(gen_0, simple_pole(gen_m2)));
2703 91 : Ldata = mkvecn(7, tag(M, t_LFUN_QF), dual,
2704 : mkvec2(gen_0, gen_1), k, d, eno, poles);
2705 91 : return gerepilecopy(ltop, Ldata);
2706 : }
2707 :
2708 : /********************************************************************/
2709 : /** Artin L function, based on a GP script by Charlotte Euvrard **/
2710 : /********************************************************************/
2711 :
2712 : static GEN
2713 119 : artin_charfromgens(GEN G, GEN M)
2714 : {
2715 119 : GEN R, V, ord = gal_get_orders(G), grp = gal_get_group(G);
2716 119 : long i, j, k, n = lg(ord)-1, m = lg(grp)-1;
2717 :
2718 119 : if (lg(M)-1 != n) pari_err_DIM("lfunartin");
2719 119 : R = cgetg(m+1, t_VEC);
2720 119 : gel(R, 1) = matid(lg(gel(M, 1))-1);
2721 357 : for (i = 1, k = 1; i <= n; ++i)
2722 : {
2723 238 : long c = k*(ord[i] - 1);
2724 238 : gel(R, ++k) = gel(M, i);
2725 1043 : for (j = 2; j <= c; ++j) gel(R, ++k) = gmul(gel(R,j), gel(M,i));
2726 : }
2727 119 : V = cgetg(m+1, t_VEC);
2728 1281 : for (i = 1; i <= m; i++) gel(V, gel(grp,i)[1]) = gtrace(gel(R,i));
2729 119 : return V;
2730 : }
2731 :
2732 : /* TODO move somewhere else? */
2733 : GEN
2734 280 : galois_get_conj(GEN G)
2735 : {
2736 280 : GEN grp = gal_get_group(G);
2737 280 : long i, k, r = lg(grp)-1;
2738 280 : GEN b = zero_F2v(r);
2739 959 : for (k = 2; k <= r; ++k)
2740 : {
2741 959 : GEN g = gel(grp,k);
2742 959 : if (!F2v_coeff(b,g[1]) && g[g[1]]==1)
2743 : {
2744 392 : pari_sp av = avma;
2745 392 : GEN F = galoisfixedfield(G, g, 1, -1);
2746 392 : if (ZX_sturmpart(F, NULL) > 0) { set_avma(av); return g; }
2747 1456 : for (i = 1; i<=r; i++)
2748 : {
2749 1344 : GEN h = gel(grp, i);
2750 1344 : long t = h[1];
2751 5264 : while (h[t]!=1) t = h[t];
2752 1344 : F2v_set(b, h[g[t]]);
2753 : }
2754 112 : set_avma(av);
2755 : }
2756 : }
2757 0 : pari_err_BUG("galois_get_conj");
2758 : return NULL;/*LCOV_EXCL_LINE*/
2759 : }
2760 :
2761 7700 : static GEN cyclotoi(GEN v) { return simplify_shallow(lift_shallow(v)); }
2762 2891 : static long cyclotos(GEN v) { return gtos(cyclotoi(v)); }
2763 2821 : static long char_dim(GEN ch) { return cyclotos(gel(ch,1)); }
2764 :
2765 : static GEN
2766 1379 : artin_gamma(GEN N, GEN G, GEN ch)
2767 : {
2768 1379 : long a, t, d = char_dim(ch);
2769 1379 : if (nf_get_r2(N) == 0) return vec01(d, 0);
2770 70 : a = galois_get_conj(G)[1];
2771 70 : t = cyclotos(gel(ch,a));
2772 70 : return vec01((d+t) / 2, (d-t) / 2);
2773 : }
2774 :
2775 : static long
2776 3213 : artin_dim(GEN ind, GEN ch)
2777 : {
2778 3213 : long n = lg(ch)-1;
2779 3213 : GEN elts = group_elts(ind, n);
2780 3213 : long i, d = lg(elts)-1;
2781 3213 : GEN s = gen_0;
2782 18123 : for(i=1; i<=d; i++)
2783 14910 : s = gadd(s, gel(ch, gel(elts,i)[1]));
2784 3213 : return gtos(gdivgu(cyclotoi(s), d));
2785 : }
2786 :
2787 : static GEN
2788 623 : artin_ind(GEN elts, GEN ch, GEN p)
2789 : {
2790 623 : long i, d = lg(elts)-1;
2791 623 : GEN s = gen_0;
2792 2149 : for(i=1; i<=d; i++)
2793 1526 : s = gadd(s, gel(ch, gmul(gel(elts,i),p)[1]));
2794 623 : return gdivgu(s, d);
2795 : }
2796 :
2797 : static GEN
2798 2282 : artin_ram(GEN nf, GEN gal, GEN aut, GEN pr, GEN ramg, GEN ch, long d)
2799 : {
2800 2282 : pari_sp av = avma;
2801 : long i, v, n;
2802 : GEN p, q, V, elts;
2803 2282 : if (d==0) return pol_1(0);
2804 616 : n = degpol(gal_get_pol(gal));
2805 616 : q = p = idealramfrobenius_aut(nf, gal, pr, ramg, aut);
2806 616 : elts = group_elts(gel(ramg,2), n);
2807 616 : v = fetch_var_higher();
2808 616 : V = cgetg(d+2, t_POL);
2809 616 : V[1] = evalsigne(1)|evalvarn(v);
2810 1239 : for(i=1; i<=d; i++)
2811 : {
2812 623 : gel(V,i+1) = artin_ind(elts, ch, q);
2813 623 : q = gmul(q, p);
2814 : }
2815 616 : delete_var();
2816 616 : V = RgXn_expint(RgX_neg(V),d+1);
2817 616 : setvarn(V,0); return gerepileupto(av, ginv(V));
2818 : }
2819 :
2820 : /* N true nf; [Artin conductor, vec of [p, Lp]] */
2821 : static GEN
2822 1379 : artin_badprimes(GEN N, GEN G, GEN aut, GEN ch)
2823 : {
2824 1379 : pari_sp av = avma;
2825 1379 : long i, d = char_dim(ch);
2826 1379 : GEN P = gel(absZ_factor(nf_get_disc(N)), 1);
2827 1379 : long lP = lg(P);
2828 1379 : GEN B = cgetg(lP, t_VEC), C = cgetg(lP, t_VEC);
2829 :
2830 3661 : for (i = 1; i < lP; ++i)
2831 : {
2832 2282 : GEN p = gel(P, i), pr = idealprimedec_galois(N, p);
2833 2282 : GEN J = idealramgroups_aut(N, G, pr, aut);
2834 2282 : GEN G0 = gel(J,2); /* inertia group */
2835 2282 : long lJ = lg(J);
2836 2282 : long sdec = artin_dim(G0, ch);
2837 2282 : long ndec = group_order(G0);
2838 2282 : long j, v = ndec * (d - sdec);
2839 3213 : for (j = 3; j < lJ; ++j)
2840 : {
2841 931 : GEN Jj = gel(J, j);
2842 931 : long s = artin_dim(Jj, ch);
2843 931 : v += group_order(Jj) * (d - s);
2844 : }
2845 2282 : gel(C, i) = powiu(p, v/ndec);
2846 2282 : gel(B, i) = mkvec2(p, artin_ram(N, G, aut, pr, J, ch, sdec));
2847 : }
2848 1379 : return gerepilecopy(av, mkvec2(ZV_prod(C), B));
2849 : }
2850 :
2851 : /* p does not divide nf.index */
2852 : static GEN
2853 52401 : idealfrobenius_easy(GEN nf, GEN gal, GEN aut, GEN T, GEN p)
2854 : {
2855 52401 : long i, l = lg(aut), f = degpol(T);
2856 52400 : GEN D, Dzk, DzkT, DXp, grp = gal_get_group(gal);
2857 52401 : pari_sp av = avma;
2858 52401 : if (f==1) return gel(grp,1);
2859 50595 : Dzk = nf_get_zkprimpart(nf);
2860 50595 : D = modii(nf_get_zkden(nf), p);
2861 50592 : DzkT = RgV_to_RgM(FqV_red(Dzk, T, p), f);
2862 50595 : DXp = RgX_to_RgC(FpX_Frobenius(T, p), f);
2863 50596 : if (!equali1(D)) DXp = FpC_Fp_mul(DXp, D, p);
2864 332160 : for(i=1; i < l; i++)
2865 : {
2866 332151 : GEN g = gel(grp,i);
2867 332151 : if (perm_orderu(g) == (ulong)f)
2868 : {
2869 170269 : GEN A = FpM_FpC_mul(DzkT, gel(aut,g[1]), p);
2870 170257 : if (ZV_equal(A, DXp)) {set_avma(av); return g; }
2871 : }
2872 : }
2873 : return NULL; /* LCOV_EXCL_LINE */
2874 : }
2875 : /* true nf; p divides nf.index, pr/p unramified */
2876 : static GEN
2877 1596 : idealfrobenius_hard(GEN nf, GEN gal, GEN aut, GEN pr)
2878 : {
2879 1596 : long i, l = lg(aut), f = pr_get_f(pr);
2880 1596 : GEN modpr, p, T, X, Xp, pi, grp = gal_get_group(gal);
2881 1596 : pari_sp av = avma;
2882 1596 : if (f==1) return gel(grp,1);
2883 1344 : pi = pr_get_gen(pr);
2884 1344 : modpr = zkmodprinit(nf, pr);
2885 1344 : p = modpr_get_p(modpr);
2886 1344 : T = modpr_get_T(modpr);
2887 1344 : X = modpr_genFq(modpr);
2888 1344 : Xp = FpX_Frobenius(T, p);
2889 9188 : for (i = 1; i < l; i++)
2890 : {
2891 9188 : GEN g = gel(grp,i);
2892 9188 : if (perm_orderu(g) == (ulong)f)
2893 : {
2894 4372 : GEN S = gel(aut,g[1]);
2895 4372 : GEN A = nf_to_Fq(nf, zk_galoisapplymod(nf,X,S,p), modpr);
2896 : /* sigma(X) = X^p (mod pr) and sigma(pi) in pr */
2897 5828 : if (ZX_equal(A, Xp) && (f == nf_get_degree(nf) ||
2898 2800 : ZC_prdvd(zk_galoisapplymod(nf,pi,S,p),pr))) { set_avma(av); return g; }
2899 : }
2900 : }
2901 : return NULL; /* LCOV_EXCL_LINE */
2902 : }
2903 :
2904 : /* true nf */
2905 : static GEN
2906 53993 : dirartin(GEN nf, GEN G, GEN V, GEN aut, GEN p, long n)
2907 : {
2908 53993 : pari_sp av = avma;
2909 : GEN pr, frob;
2910 : /* pick one maximal ideal in the conjugacy class above p */
2911 53993 : GEN T = nf_get_pol(nf);
2912 53993 : if (!dvdii(nf_get_index(nf), p))
2913 : { /* simple case */
2914 52396 : GEN F = FpX_factor(T, p), P = gmael(F,1,1);
2915 52402 : frob = idealfrobenius_easy(nf, G, aut, P, p);
2916 : }
2917 : else
2918 : {
2919 1596 : pr = idealprimedec_galois(nf,p);
2920 1596 : frob = idealfrobenius_hard(nf, G, aut, pr);
2921 : }
2922 53997 : set_avma(av); return n ? RgXn_inv(gel(V, frob[1]), n): gel(V, frob[1]);
2923 : }
2924 :
2925 : GEN
2926 15666 : dirartin_worker(GEN P, ulong X, GEN nf, GEN G, GEN V, GEN aut)
2927 : {
2928 15666 : pari_sp av = avma;
2929 15666 : long i, l = lg(P);
2930 15666 : GEN W = cgetg(l, t_VEC);
2931 69629 : for(i = 1; i < l; i++)
2932 : {
2933 53963 : ulong p = uel(P,i);
2934 53963 : long d = ulogint(X, p) + 1; /* minimal d such that p^d > X */
2935 53965 : gel(W,i) = dirartin(nf, G, V, aut, utoi(uel(P,i)), d);
2936 : }
2937 15666 : return gerepilecopy(av, mkvec2(P,W));
2938 : }
2939 :
2940 : static GEN
2941 2947 : vecan_artin(GEN an, long L, long prec)
2942 : {
2943 2947 : GEN A, Sbad = gel(an,5);
2944 2947 : long n = itos(gel(an,6)), isreal = lg(an)<8 ? 0: !itos(gel(an,7));
2945 2947 : GEN worker = snm_closure(is_entry("_dirartin_worker"), vecslice(an,1,4));
2946 2947 : A = lift_shallow(pardireuler(worker, gen_2, stoi(L), NULL, Sbad));
2947 2947 : A = RgXV_RgV_eval(A, grootsof1(n, prec));
2948 2947 : if (isreal) A = real_i(A);
2949 2947 : return A;
2950 : }
2951 :
2952 : static GEN
2953 28 : eulerf_artin(GEN an, GEN p, long prec)
2954 : {
2955 28 : GEN nf = gel(an,1), G = gel(an,2), V = gel(an,3), aut = gel(an,4);
2956 28 : GEN Sbad = gel(an,5);
2957 28 : long n = itos(gel(an,6)), isreal = lg(an)<8 ? 0: !itos(gel(an,7));
2958 28 : GEN f = eulerf_bad(Sbad, p);
2959 28 : if (!f) f = mkrfrac(gen_1,dirartin(nf, G, V, aut, p, 0));
2960 28 : f = gsubst(liftpol(f),1, rootsof1u_cx(n, prec));
2961 28 : if (isreal) f = real_i(f);
2962 28 : return f;
2963 : }
2964 :
2965 : static GEN
2966 2856 : char_expand(GEN conj, GEN ch)
2967 : {
2968 2856 : long i, l = lg(conj);
2969 2856 : GEN V = cgetg(l, t_COL);
2970 31409 : for (i=1; i<l; i++) gel(V,i) = gel(ch, conj[i]);
2971 2856 : return V;
2972 : }
2973 :
2974 : static GEN
2975 1596 : handle_zeta(long n, GEN ch, long *m)
2976 : {
2977 : GEN c;
2978 1596 : long t, i, l = lg(ch);
2979 1596 : GEN dim = cyclotoi(vecsum(ch));
2980 1596 : if (typ(dim) != t_INT)
2981 0 : pari_err_DOMAIN("lfunartin","chi","is not a", strtoGENstr("character"), ch);
2982 1596 : t = itos(dim);
2983 1596 : if (t < 0 || t % n)
2984 0 : pari_err_DOMAIN("lfunartin","chi","is not a", strtoGENstr("character"), ch);
2985 1596 : if (t == 0) { *m = 0; return ch; }
2986 224 : *m = t / n;
2987 224 : c = cgetg(l, t_COL);
2988 2065 : for (i=1; i<l; i++)
2989 1841 : gel(c,i) = gsubgs(gel(ch,i), *m);
2990 224 : return c;
2991 : }
2992 :
2993 : static int
2994 6496 : cyclo_is_real(GEN v, GEN ix)
2995 : {
2996 6496 : pari_sp av = avma;
2997 6496 : GEN w = poleval(lift_shallow(v), ix);
2998 6496 : return gc_bool(av, gequal(w, v));
2999 : }
3000 :
3001 : static int
3002 1379 : char_is_real(GEN ch, GEN mod)
3003 : {
3004 1379 : long i, l = lg(ch);
3005 1379 : GEN ix = QXQ_inv(pol_x(varn(mod)), mod);
3006 7014 : for (i=1; i<l; i++)
3007 6496 : if (!cyclo_is_real(gel(ch,i), ix)) return 0;
3008 518 : return 1;
3009 : }
3010 :
3011 : GEN
3012 1610 : lfunartin(GEN nf, GEN gal, GEN ch, long o, long bitprec)
3013 : {
3014 1610 : pari_sp av = avma;
3015 1610 : GEN bc, V, aut, mod, Ldata = NULL, chx, cc, conj, repr;
3016 : long tmult, var;
3017 1610 : nf = checknf(nf);
3018 1610 : checkgal(gal);
3019 1610 : var = gvar(ch);
3020 1610 : if (var == 0) pari_err_PRIORITY("lfunartin",ch,"=",0);
3021 1610 : if (var < 0) var = 1;
3022 1610 : if (!is_vec_t(typ(ch))) pari_err_TYPE("lfunartin", ch);
3023 1610 : cc = group_to_cc(gal);
3024 1610 : conj = gel(cc,2);
3025 1610 : repr = gel(cc,3);
3026 1610 : mod = mkpolmod(gen_1, polcyclo(o, var));
3027 1610 : if (lg(ch)>1 && typ(gel(ch,1))==t_MAT)
3028 119 : chx = artin_charfromgens(gal, gmul(ch,mod));
3029 : else
3030 : {
3031 1491 : if (lg(repr) != lg(ch)) pari_err_DIM("lfunartin");
3032 1477 : chx = char_expand(conj, gmul(ch,mod));
3033 : }
3034 1596 : chx = handle_zeta(nf_get_degree(nf), chx, &tmult);
3035 1596 : ch = shallowextract(chx, repr);
3036 1596 : if (!gequal0(chx))
3037 : {
3038 1379 : GEN real = char_is_real(chx, gel(mod,1))? gen_0: gen_1;
3039 1379 : aut = nfgaloispermtobasis(nf, gal);
3040 1379 : V = gmul(char_expand(conj, galoischarpoly(gal, ch, o)), mod);
3041 1379 : bc = artin_badprimes(nf, gal, aut, chx);
3042 2758 : Ldata = mkvecn(6,
3043 1379 : tag(mkcoln(7, nf, gal, V, aut, gel(bc, 2), stoi(o), real), t_LFUN_ARTIN),
3044 1379 : real, artin_gamma(nf, gal, chx), gen_1, gel(bc,1), gen_0);
3045 : }
3046 1596 : if (tmult==0 && Ldata==NULL) pari_err_TYPE("lfunartin",ch);
3047 1596 : if (tmult)
3048 : {
3049 : long i;
3050 224 : if (Ldata==NULL) { Ldata = lfunzeta(); tmult--; }
3051 231 : for(i=1; i<=tmult; i++)
3052 7 : Ldata = lfunmul(Ldata, gen_1, bitprec);
3053 : }
3054 1596 : return gerepilecopy(av, Ldata);
3055 : }
3056 :
3057 : /* true nf */
3058 : static GEN
3059 21 : lfunzetakinit_artin(GEN nf, GEN gal, GEN dom, long der, long bitprec)
3060 : {
3061 21 : GEN F, E, M, domain, To = galoischartable(gal), T = gel(To, 1);
3062 21 : long i, o = itos(gel(To, 2)), l = lg(T);
3063 21 : F = cgetg(l, t_VEC);
3064 21 : E = cgetg(l, t_VECSMALL);
3065 84 : for (i = 1; i < l; ++i)
3066 : {
3067 63 : GEN L = lfunartin(nf, gal, gel(T,i), o, bitprec);
3068 63 : gel(F, i) = lfuninit(L, dom, der, bitprec);
3069 63 : E[i] = char_dim(gel(T,i));
3070 : }
3071 21 : domain = mkvec2(dom, mkvecsmall2(der, bitprec));
3072 21 : M = mkvec3(F, E, zero_zv(l-1));
3073 21 : return lfuninit_make(t_LDESC_PRODUCT, lfunzetak_i(nf), M, domain);
3074 : }
3075 :
3076 : /********************************************************************/
3077 : /** High-level Constructors **/
3078 : /********************************************************************/
3079 : enum { t_LFUNMISC_POL, t_LFUNMISC_CHIQUAD, t_LFUNMISC_CHICONREY,
3080 : t_LFUNMISC_CHIGEN, t_LFUNMISC_ELLINIT, t_LFUNMISC_ETAQUO,
3081 : t_LFUNMISC_GCHAR, t_LFUNMISC_ABELREL };
3082 : static long
3083 17696 : lfundatatype(GEN data)
3084 : {
3085 17696 : switch(typ(data))
3086 : {
3087 3899 : case t_INT: return t_LFUNMISC_CHIQUAD;
3088 217 : case t_INTMOD: return t_LFUNMISC_CHICONREY;
3089 651 : case t_POL: return t_LFUNMISC_POL;
3090 12929 : case t_VEC:
3091 12929 : switch(lg(data))
3092 : {
3093 4053 : case 17: return t_LFUNMISC_ELLINIT;
3094 0 : case 10: return t_LFUNMISC_POL;
3095 8876 : case 3:
3096 8876 : if (typ(gel(data,1)) != t_VEC) break;
3097 8876 : return is_gchar_group(gel(data,1)) ? t_LFUNMISC_GCHAR
3098 8876 : : typ(gel(data,2))==t_MAT ? t_LFUNMISC_ABELREL
3099 : : t_LFUNMISC_CHIGEN;
3100 : }
3101 0 : break;
3102 : }
3103 0 : return -1;
3104 : }
3105 : static GEN
3106 122114 : lfunmisc_to_ldata_i(GEN ldata, long shallow)
3107 : {
3108 : GEN x;
3109 122114 : if (is_linit(ldata)) ldata = linit_get_ldata(ldata);
3110 122114 : if (is_ldata(ldata) && is_tagged(ldata))
3111 : {
3112 104250 : if (!shallow) ldata = gcopy(ldata);
3113 104250 : checkldata(ldata); return ldata;
3114 : }
3115 17864 : x = checknf_i(ldata); if (x) return lfunzetak(x);
3116 17696 : switch (lfundatatype(ldata))
3117 : {
3118 651 : case t_LFUNMISC_POL: return lfunzetak(ldata);
3119 3899 : case t_LFUNMISC_CHIQUAD: return lfunchiquad(ldata);
3120 217 : case t_LFUNMISC_CHICONREY:
3121 : {
3122 217 : GEN G = znstar0(gel(ldata,1), 1);
3123 217 : return lfunchiZ(G, gel(ldata,2));
3124 : }
3125 8463 : case t_LFUNMISC_CHIGEN:
3126 : {
3127 8463 : GEN G = gel(ldata,1), chi = gel(ldata,2);
3128 8463 : switch(nftyp(G))
3129 : {
3130 8204 : case typ_BIDZ: return lfunchiZ(G, chi);
3131 252 : case typ_BNR: return lfunchigen(G, chi);
3132 : }
3133 : }
3134 7 : break;
3135 392 : case t_LFUNMISC_GCHAR: return lfungchar(gel(ldata,1), gel(ldata,2));
3136 21 : case t_LFUNMISC_ABELREL:
3137 21 : return lfunabelrel(gel(ldata,1), gel(ldata,2),
3138 21 : lg(ldata)==3? NULL: gel(ldata,3));
3139 4053 : case t_LFUNMISC_ELLINIT: return lfunell(ldata);
3140 : }
3141 7 : if (shallow != 2) pari_err_TYPE("lfunmisc_to_ldata",ldata);
3142 0 : return NULL;
3143 : }
3144 :
3145 : GEN
3146 3514 : lfunmisc_to_ldata(GEN ldata)
3147 3514 : { return lfunmisc_to_ldata_i(ldata, 0); }
3148 :
3149 : GEN
3150 97803 : lfunmisc_to_ldata_shallow(GEN ldata)
3151 97803 : { return lfunmisc_to_ldata_i(ldata, 1); }
3152 :
3153 : GEN
3154 20797 : lfunmisc_to_ldata_shallow_i(GEN ldata)
3155 20797 : { return lfunmisc_to_ldata_i(ldata, 2); }
3156 :
3157 : /********************************************************************/
3158 : /** High-level an expansion **/
3159 : /********************************************************************/
3160 : /* van is the output of ldata_get_an: return a_1,...a_L at precision prec */
3161 : GEN
3162 31878 : ldata_vecan(GEN van, long L, long prec)
3163 : {
3164 31878 : GEN an = gel(van, 2);
3165 31878 : long t = mael(van,1,1);
3166 : pari_timer ti;
3167 31878 : if (DEBUGLEVEL >= 1)
3168 0 : err_printf("Lfun: computing %ld coeffs, prec %ld, type %ld\n", L, prec, t);
3169 31878 : if (DEBUGLEVEL >= 2) timer_start(&ti);
3170 31878 : if (L < 0) L = 0;
3171 31878 : switch (t)
3172 : {
3173 : long n;
3174 1932 : case t_LFUN_GENERIC:
3175 1932 : an = vecan_closure(an, L, prec);
3176 1862 : n = lg(an)-1;
3177 1862 : if (n < L)
3178 : {
3179 14 : pari_warn(warner, "#an = %ld < %ld, results may be imprecise", n, L);
3180 14 : an = shallowconcat(an, zerovec(L-n));
3181 : }
3182 1862 : break;
3183 0 : case t_LFUN_CLOSURE0:
3184 : pari_err_BUG("ldata_vecan: please call ldata_newprec");/*LCOV_EXCL_LINE*/
3185 2898 : case t_LFUN_ZETA: an = const_vecsmall(L, 1); break;
3186 2100 : case t_LFUN_NF: an = dirzetak(an, stoi(L)); break;
3187 6258 : case t_LFUN_ELL:
3188 6258 : an = (ell_get_type(an) == t_ELL_Q) ? ellanQ_zv(an, L): ellan(an, L);
3189 6258 : break;
3190 3654 : case t_LFUN_KRONECKER: an = vecan_Kronecker(an, L); break;
3191 168 : case t_LFUN_ABELREL: an = vecan_abelrel(an, L); break;
3192 3647 : case t_LFUN_CHIZ: an = vecan_chiZ(an, L, prec); break;
3193 1344 : case t_LFUN_CHIGEN: an = vecan_chigen(an, L, prec); break;
3194 693 : case t_LFUN_HECKE: an = vecan_gchar(an, L, prec); break;
3195 2947 : case t_LFUN_ARTIN: an = vecan_artin(an, L, prec); break;
3196 147 : case t_LFUN_ETA: an = vecan_eta(an, L); break;
3197 399 : case t_LFUN_QF: an = vecan_qf(an, L); break;
3198 630 : case t_LFUN_DIV: an = vecan_div(an, L, prec); break;
3199 308 : case t_LFUN_MUL: an = vecan_mul(an, L, prec); break;
3200 126 : case t_LFUN_CONJ: an = vecan_conj(an, L, prec); break;
3201 343 : case t_LFUN_SYMPOW_ELL: an = vecan_ellsympow(an, L); break;
3202 196 : case t_LFUN_GENUS2: an = vecan_genus2(an, L); break;
3203 168 : case t_LFUN_HGM:
3204 168 : an = hgmcoefs(gel(an,1), gel(an,2), L); break;
3205 406 : case t_LFUN_MFCLOS:
3206 : {
3207 406 : GEN F = gel(an,1), E = gel(an,2), c = gel(an,3);
3208 406 : an = mfcoefs(F,L,1) + 1; /* skip a_0 */
3209 406 : an[0] = evaltyp(t_VEC)|_evallg(L+1);
3210 406 : an = mfvecembed(E, an);
3211 406 : if (!isint1(c)) an = RgV_Rg_mul(an,c);
3212 406 : break;
3213 : }
3214 2877 : case t_LFUN_TWIST: an = vecan_twist(an, L, prec); break;
3215 637 : case t_LFUN_SHIFT: an = vecan_shift(an, L, prec); break;
3216 0 : default: pari_err_TYPE("ldata_vecan", van);
3217 : }
3218 31808 : if (DEBUGLEVEL >= 2) timer_printf(&ti, "ldata_vecan");
3219 31808 : return an;
3220 : }
3221 :
3222 : /* shallow */
3223 : GEN
3224 32613 : ldata_newprec(GEN ldata, long prec)
3225 : {
3226 32613 : GEN van = ldata_get_an(ldata), an = gel(van, 2);
3227 32613 : long t = mael(van,1,1);
3228 32613 : switch (t)
3229 : {
3230 154 : case t_LFUN_CLOSURE0: return closure2ldata(an, prec);
3231 994 : case t_LFUN_HECKE:
3232 : {
3233 994 : GEN gc = gel(an, 1), chiw = gel(an, 2);
3234 994 : gc = gcharnewprec(gc, prec);
3235 994 : return gchari_lfun(gc, chiw, gen_0); /* chi in internal format */
3236 : }
3237 329 : case t_LFUN_QF:
3238 : {
3239 329 : GEN eno = ldata_get_rootno(ldata);
3240 329 : if (typ(eno)==t_REAL && realprec(eno) < prec) return lfunqf(an, prec);
3241 273 : break;
3242 : }
3243 : }
3244 31409 : return ldata;
3245 : }
3246 :
3247 : GEN
3248 854 : ldata_eulerf(GEN van, GEN p, long prec)
3249 : {
3250 854 : GEN an = gel(van, 2), f = gen_0;
3251 854 : long t = mael(van,1,1);
3252 854 : switch (t)
3253 : {
3254 70 : case t_LFUN_GENERIC:
3255 70 : f = eulerf_closure(an, p, prec); break;
3256 0 : case t_LFUN_CLOSURE0:
3257 : pari_err_BUG("ldata_vecan: please call ldata_newprec");/*LCOV_EXCL_LINE*/
3258 70 : case t_LFUN_ZETA: f = mkrfrac(gen_1,deg1pol(gen_m1, gen_1,0)); break;
3259 168 : case t_LFUN_NF: f = eulerf_zetak(an, p); break;
3260 70 : case t_LFUN_ELL: f = elleulerf(an, p); break;
3261 70 : case t_LFUN_KRONECKER:
3262 70 : f = mkrfrac(gen_1, deg1pol_shallow(stoi(-kronecker(an, p)), gen_1, 0)); break;
3263 42 : case t_LFUN_ABELREL: f = eulerf_abelrel(an, p); break;
3264 42 : case t_LFUN_CHIZ: f = eulerf_chiZ(an, p, prec); break;
3265 28 : case t_LFUN_CHIGEN: f = eulerf_chigen(an, p, prec); break;
3266 14 : case t_LFUN_HECKE: f = eulerf_gchar(an, p, prec); break;
3267 28 : case t_LFUN_ARTIN: f = eulerf_artin(an, p, prec); break;
3268 14 : case t_LFUN_DIV: f = eulerf_div(an, p, prec); break;
3269 28 : case t_LFUN_MUL: f = eulerf_mul(an, p, prec); break;
3270 0 : case t_LFUN_CONJ: f = eulerf_conj(an, p, prec); break;
3271 28 : case t_LFUN_SYMPOW_ELL: f = eulerf_ellsympow(an, p); break;
3272 0 : case t_LFUN_GENUS2: f = eulerf_genus2(an, p); break;
3273 14 : case t_LFUN_TWIST: f = eulerf_twist(an, p, prec); break;
3274 42 : case t_LFUN_SHIFT: f = eulerf_shift(an, p, prec); break;
3275 84 : case t_LFUN_HGM: f = eulerf_hgm(an, p); break;
3276 42 : default: f = NULL; break;
3277 : }
3278 854 : if (!f) pari_err_DOMAIN("lfuneuler", "L", "Euler product", strtoGENstr("unknown"), an);
3279 742 : return f;
3280 : }
3281 :
3282 : GEN
3283 707 : lfuneuler(GEN ldata, GEN p, long prec)
3284 : {
3285 707 : pari_sp av = avma;
3286 707 : if (typ(p)!=t_INT || signe(p)<=0) pari_err_TYPE("lfuneuler", p);
3287 700 : ldata = ldata_newprec(lfunmisc_to_ldata_shallow(ldata), prec);
3288 700 : return gerepilecopy(av, ldata_eulerf(ldata_get_an(ldata), p, prec));
3289 : }
|
