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 **/
18 : /** **/
19 : /********************************************************************/
20 :
21 : #include "pari.h"
22 : #include "paripriv.h"
23 :
24 : #define DEBUGLEVEL DEBUGLEVEL_lfun
25 :
26 : /*******************************************************************/
27 : /* Accessors */
28 : /*******************************************************************/
29 :
30 : static GEN
31 12913 : mysercoeff(GEN x, long n)
32 : {
33 12913 : long N = n - valser(x);
34 12913 : return (N < 0)? gen_0: gel(x, N+2);
35 : }
36 :
37 : long
38 78290 : ldata_get_type(GEN ldata) { return mael3(ldata, 1, 1, 1); }
39 :
40 : GEN
41 77106 : ldata_get_an(GEN ldata) { return gel(ldata, 1); }
42 :
43 : GEN
44 61495 : ldata_get_dual(GEN ldata) { return gel(ldata, 2); }
45 :
46 : long
47 2834 : ldata_isreal(GEN ldata) { return isintzero(gel(ldata, 2)); }
48 :
49 : GEN
50 357148 : ldata_get_gammavec(GEN ldata) { return gel(ldata, 3); }
51 :
52 : long
53 26157 : ldata_get_degree(GEN ldata) { return lg(gel(ldata, 3))-1; }
54 :
55 : GEN
56 169561 : ldata_get_k(GEN ldata)
57 : {
58 169561 : GEN w = gel(ldata,4);
59 169561 : if (typ(w) == t_VEC) w = gel(w,1);
60 169561 : return w;
61 : }
62 :
63 : /* a_n = O(n^{k1 + epsilon}) */
64 : GEN
65 224 : ldata_get_k1(GEN ldata)
66 : {
67 224 : GEN w = gel(ldata,4);
68 224 : if (typ(w) == t_VEC) return gel(w,2);
69 : /* by default, assume that k1 = k-1 and even (k-1)/2 for entire functions */
70 224 : w = gaddgs(w,-1);
71 224 : return ldata_get_residue(ldata)? w: gmul2n(w, -1);
72 : }
73 :
74 : /* a_n = O(n^{k1 + epsilon}) */
75 : static double
76 93172 : ldata_get_k1_dbl(GEN ldata)
77 : {
78 93172 : GEN w = gel(ldata,4);
79 : double k;
80 93172 : if (typ(w) == t_VEC) return gtodouble(gel(w,2));
81 : /* by default, assume that k1 = k-1 and even (k-1)/2 for entire functions */
82 91527 : k = gtodouble(w);
83 91527 : return ldata_get_residue(ldata)? k-1: (k-1)/2.;
84 : }
85 :
86 : GEN
87 279040 : ldata_get_conductor(GEN ldata) { return gel(ldata, 5); }
88 :
89 : GEN
90 108318 : ldata_get_rootno(GEN ldata) { return gel(ldata, 6); }
91 :
92 : GEN
93 179956 : ldata_get_residue(GEN ldata) { return lg(ldata) == 7 ? NULL: gel(ldata, 7); }
94 :
95 : long
96 149139 : linit_get_type(GEN linit) { return mael(linit, 1, 1); }
97 :
98 : GEN
99 200312 : linit_get_ldata(GEN linit) { return gel(linit, 2); }
100 :
101 : GEN
102 252802 : linit_get_tech(GEN linit) { return gel(linit, 3); }
103 :
104 : long
105 305828 : is_linit(GEN data)
106 : {
107 186165 : return lg(data) == 4 && typ(data) == t_VEC
108 491993 : && typ(gel(data, 1)) == t_VECSMALL;
109 : }
110 :
111 : GEN
112 32513 : lfun_get_step(GEN tech) { return gmael(tech, 2, 1);}
113 :
114 : GEN
115 32513 : lfun_get_pol(GEN tech) { return gmael(tech, 2, 2);}
116 :
117 : GEN
118 5431 : lfun_get_Residue(GEN tech) { return gmael(tech, 2, 3);}
119 :
120 : GEN
121 50628 : lfun_get_k2(GEN tech) { return gmael(tech, 3, 1);}
122 :
123 : GEN
124 19382 : lfun_get_w2(GEN tech) { return gmael(tech, 3, 2);}
125 :
126 : GEN
127 19382 : lfun_get_expot(GEN tech) { return gmael(tech, 3, 3);}
128 :
129 : GEN
130 10402 : lfun_get_factgammavec(GEN tech) { return gmael(tech, 3, 4); }
131 :
132 : /* Handle complex Vga whose sum is real */
133 : static GEN
134 105527 : sumVga(GEN Vga) { return real_i(vecsum(Vga)); }
135 : /* sum_i max (Im v[i],0) */
136 : static double
137 27617 : sumVgaimpos(GEN v)
138 : {
139 27617 : double d = 0.;
140 27617 : long i, l = lg(v);
141 76845 : for (i = 1; i < l; i++)
142 : {
143 49228 : GEN c = imag_i(gel(v,i));
144 49228 : if (gsigne(c) > 0) d += gtodouble(c);
145 : }
146 27617 : return d;
147 : }
148 :
149 : static long
150 44432 : vgaell(GEN Vga)
151 : {
152 44432 : if (lg(Vga) == 3)
153 31209 : { GEN c = gsub(gel(Vga,1), gel(Vga,2)); return gequal1(c) || gequalm1(c); }
154 13223 : return 0;
155 : }
156 : int
157 88851 : Vgaeasytheta(GEN Vga) { return lg(Vga)-1 == 1 || vgaell(Vga); }
158 : /* return b(n) := a(n) * n^c, when Vgaeasytheta(Vga) is set */
159 : static GEN
160 18725 : antwist(GEN an, GEN Vga, long prec)
161 : {
162 : long l, i;
163 18725 : GEN b, c = vecmin(Vga);
164 18725 : if (gequal0(c)) return an;
165 4466 : l = lg(an); b = cgetg(l, t_VEC);
166 4466 : if (gequal1(c))
167 : {
168 3626 : if (typ(an) == t_VECSMALL)
169 17647 : for (i = 1; i < l; i++) gel(b,i) = mulss(an[i], i);
170 : else
171 41356 : for (i = 1; i < l; i++) gel(b,i) = gmulgu(gel(an,i), i);
172 : }
173 : else
174 : {
175 840 : GEN v = vecpowug(l-1, c, prec);
176 840 : if (typ(an) == t_VECSMALL)
177 0 : for (i = 1; i < l; i++) gel(b,i) = gmulsg(an[i], gel(v,i));
178 : else
179 34573 : for (i = 1; i < l; i++) gel(b,i) = gmul(gel(an,i), gel(v,i));
180 : }
181 4466 : return b;
182 : }
183 :
184 : static GEN
185 10150 : theta_dual(GEN theta, GEN bn)
186 : {
187 10150 : if (typ(bn)==t_INT) return NULL;
188 : else
189 : {
190 77 : GEN thetad = shallowcopy(theta), ldata = linit_get_ldata(theta);
191 77 : GEN Vga = ldata_get_gammavec(ldata);
192 77 : GEN tech = shallowcopy(linit_get_tech(theta));
193 77 : GEN an = theta_get_an(tech);
194 77 : long prec = nbits2prec(theta_get_bitprec(tech));
195 77 : GEN vb = ldata_vecan(bn, lg(an)-1, prec);
196 77 : if (!theta_get_m(tech) && Vgaeasytheta(Vga)) vb = antwist(vb, Vga, prec);
197 77 : gel(tech,1) = vb;
198 77 : gel(thetad,3) = tech; return thetad;
199 : }
200 : }
201 :
202 : static GEN
203 85420 : domain_get_dom(GEN domain) { return gel(domain,1); }
204 : static long
205 25708 : domain_get_der(GEN domain) { return mael2(domain, 2, 1); }
206 : static long
207 37517 : domain_get_bitprec(GEN domain) { return mael2(domain, 2, 2); }
208 : GEN
209 86001 : lfun_get_domain(GEN tech) { return gel(tech,1); }
210 : long
211 91 : lfun_get_bitprec(GEN tech){ return domain_get_bitprec(lfun_get_domain(tech)); }
212 : GEN
213 60125 : lfun_get_dom(GEN tech) { return domain_get_dom(lfun_get_domain(tech)); }
214 :
215 : GEN
216 2575 : lfunprod_get_fact(GEN tech) { return gel(tech, 2); }
217 :
218 : GEN
219 52780 : theta_get_an(GEN tdata) { return gel(tdata, 1);}
220 : GEN
221 9002 : theta_get_K(GEN tdata) { return gel(tdata, 2);}
222 : GEN
223 5957 : theta_get_R(GEN tdata) { return gel(tdata, 3);}
224 : long
225 66668 : theta_get_bitprec(GEN tdata) { return itos(gel(tdata, 4));}
226 : long
227 102641 : theta_get_m(GEN tdata) { return itos(gel(tdata, 5));}
228 : GEN
229 54754 : theta_get_tdom(GEN tdata) { return gel(tdata, 6);}
230 : GEN
231 63735 : theta_get_isqrtN(GEN tdata) { return gel(tdata, 7);}
232 :
233 : /*******************************************************************/
234 : /* Helper functions related to Gamma products */
235 : /*******************************************************************/
236 : /* x != 0 */
237 : static int
238 7147 : serisscalar(GEN x)
239 : {
240 : long i;
241 7147 : if (valser(x)) return 0;
242 9520 : for (i = lg(x)-1; i > 3; i--) if (!gequal0(gel(x,i))) return 0;
243 6895 : return 1;
244 : }
245 :
246 : /* return -itos(s) >= 0 if scalar s is (approximately) equal to a nonpositive
247 : * integer, and -1 otherwise */
248 : static long
249 22127 : isnegint(GEN s)
250 : {
251 22127 : GEN r = ground(real_i(s));
252 22127 : if (signe(r) <= 0 && gequal(s, r)) return -itos(r);
253 22001 : return -1;
254 : }
255 : /* if s = a + O(x^n), a <= 0 integer, replace by a + b*x^n + O(x^(n+1)) */
256 : static GEN
257 7168 : serextendifnegint(GEN s, GEN b, long *ext)
258 : {
259 7168 : if (!signe(s) || (serisscalar(s) && isnegint(gel(s,2)) >= 0))
260 : {
261 112 : long l = lg(s);
262 112 : GEN t = cgetg(l+1, t_SER);
263 301 : gel(t, l) = b; while (--l > 1) gel(t,l) = gel(s,l);
264 112 : if (gequal0(gel(t,2))) gel(t,2) = gen_0;
265 112 : t[1] = s[1]; s = normalizeser(t); *ext = 1;
266 : }
267 7168 : return s;
268 : }
269 :
270 : /* r/x + O(1), r != 0 */
271 : static GEN
272 5131 : serpole(GEN r)
273 : {
274 5131 : GEN s = cgetg(3, t_SER);
275 5131 : s[1] = evalsigne(1)|evalvalser(-1)|evalvarn(0);
276 5131 : gel(s,2) = r; return s;
277 : }
278 : /* a0 + a1 x + O(x^e), e >= 0 */
279 : static GEN
280 8575 : deg1ser_shallow(GEN a1, GEN a0, long v, long e)
281 8575 : { return RgX_to_ser(deg1pol_shallow(a1, a0, v), e+2); }
282 :
283 : /* pi^(-s/2) Gamma(s/2) */
284 : static GEN
285 10934 : gamma_R(GEN s, long *ext, long prec)
286 : {
287 10934 : GEN s2 = gmul2n(s, -1);
288 : long ms;
289 :
290 10934 : if (typ(s) == t_SER)
291 5068 : s2 = serextendifnegint(s2, ghalf, ext);
292 5866 : else if ((ms = isnegint(s2)) >= 0)
293 : {
294 35 : GEN r = gmul(powPis(stoi(ms),prec), gdivsg(odd(ms)? -2: 2, mpfact(ms)));
295 35 : return serpole(r);
296 : }
297 10899 : return gdiv(ggamma(s2,prec), powPis(s2,prec));
298 : }
299 : /* gamma_R(s)gamma_R(s+1) = 2 (2pi)^(-s) Gamma(s) */
300 : static GEN
301 11466 : gamma_C(GEN s, long *ext, long prec)
302 : {
303 : long ms;
304 11466 : if (typ(s) == t_SER)
305 2100 : s = serextendifnegint(s, gen_1, ext);
306 9366 : else if ((ms = isnegint(s)) >= 0)
307 : {
308 0 : GEN r = gmul(pow2Pis(stoi(ms),prec), gdivsg(odd(ms)? -2: 2, mpfact(ms)));
309 0 : return serpole(r);
310 : }
311 11466 : return gmul2n(gdiv(ggamma(s,prec), pow2Pis(s,prec)), 1);
312 : }
313 :
314 : static GEN
315 2254 : gammafrac(GEN r, long d)
316 : {
317 2254 : long i, l = labs(d) + 1, j = (d > 0)? 0: 2*d;
318 2254 : GEN T, v = cgetg(l, t_COL);
319 7385 : for (i = 1; i < l; i++, j += 2)
320 5131 : gel(v,i) = deg1pol_shallow(gen_1, gaddgs(r, j), 0);
321 2254 : T = RgV_prod(v); return d > 0? T: mkrfrac(gen_1, T);
322 : }
323 :
324 : /*
325 : GR(s)=Pi^-(s/2)*gamma(s/2);
326 : GC(s)=2*(2*Pi)^-s*gamma(s)
327 : gdirect(F,s)=prod(i=1,#F,GR(s+F[i]))
328 : gfact(F,s)=
329 : { my([R,A,B]=gammafactor(F), [a,e]=A, [b,f]=B, p=poldegree(R));
330 : subst(R,x,s) * (2*Pi)^-p * prod(i=1,#a,GR(s+a[i])^e[i])
331 : * prod(i=1,#b,GC(s+b[i])^f[i]); }
332 : */
333 : static GEN
334 22841 : gammafactor(GEN Vga)
335 : {
336 22841 : long i, r, c, l = lg(Vga);
337 22841 : GEN v, P, a, b, e, f, E, F = cgetg(l, t_VEC), R = gen_1;
338 64792 : for (i = 1; i < l; ++i)
339 : {
340 41951 : GEN a = gel(Vga,i), r = gmul2n(real_i(a), -1);
341 41951 : long q = itos(gfloor(r)); /* [Re a/2] */
342 41951 : r = gmul2n(gsubgs(r, q), 1);
343 41951 : gel(F,i) = gequal0(imag_i(a)) ? r : mkcomplex(r, imag_i(a)); /* 2{Re a/2} + I*(Im a) */
344 41951 : if (q) R = gmul(R, gammafrac(gel(F,i), q));
345 : }
346 22841 : F = vec_reduce(F, &E); l = lg(E);
347 22841 : v = cgetg(l, t_VEC);
348 58114 : for (i = 1; i < l; i++)
349 35273 : gel(v,i) = mkvec2(gsub(gel(F,i),gfloor(real_i(gel(F,i)))), stoi(E[i]));
350 22841 : gen_sort_inplace(v, (void*)cmp_universal, cmp_nodata, &P);
351 22841 : a = cgetg(l, t_VEC); e = cgetg(l, t_VECSMALL);
352 22841 : b = cgetg(l, t_VEC); f = cgetg(l, t_VECSMALL);
353 46774 : for (i = r = c = 1; i < l;)
354 23933 : if (i==l-1 || cmp_universal(gel(v,i), gel(v,i+1)))
355 12593 : { gel(a, r) = gel(F, P[i]); e[r++] = E[P[i]]; i++; }
356 : else
357 11340 : { gel(b, c) = gel(F, P[i]); f[c++] = E[P[i]]; i+=2; }
358 22841 : setlg(a, r); setlg(e, r);
359 22841 : setlg(b, c); setlg(f, c); return mkvec3(R, mkvec2(a,e), mkvec2(b,f));
360 : }
361 :
362 : static GEN
363 5054 : polgammaeval(GEN F, GEN s)
364 : {
365 5054 : GEN r = poleval(F, s);
366 5054 : if (typ(s) != t_SER && gequal0(r))
367 : { /* here typ(F) = t_POL */
368 : long e;
369 7 : for (e = 1;; e++)
370 : {
371 7 : F = RgX_deriv(F); r = poleval(F,s);
372 7 : if (!gequal0(r)) break;
373 : }
374 7 : if (e > 1) r = gdiv(r, mpfact(e));
375 7 : r = serpole(r); setvalser(r, e);
376 : }
377 5054 : return r;
378 : }
379 : static long
380 2499 : rfrac_degree(GEN R)
381 : {
382 2499 : GEN a = gel(R,1), b = gel(R,2);
383 2499 : return ((typ(a) == t_POL)? degpol(a): 0) - degpol(b);
384 : }
385 : static GEN
386 21280 : fracgammaeval(GEN F, GEN s, long prec)
387 : {
388 21280 : GEN R = gel(F,1);
389 : long d;
390 21280 : switch(typ(R))
391 : {
392 56 : case t_POL:
393 56 : d = degpol(R);
394 56 : R = polgammaeval(R, s); break;
395 2499 : case t_RFRAC:
396 2499 : d = rfrac_degree(R);
397 2499 : R = gdiv(polgammaeval(gel(R,1), s), polgammaeval(gel(R,2), s)); break;
398 18725 : default: return R;
399 : }
400 2555 : return gmul(R, powrs(Pi2n(1,prec), -d));
401 : }
402 :
403 : static GEN
404 21280 : gammafactproduct(GEN F, GEN s, long *ext, long prec)
405 : {
406 21280 : pari_sp av = avma;
407 21280 : GEN R = gel(F,2), Rw = gel(R,1), Re = gel(R,2);
408 21280 : GEN C = gel(F,3), Cw = gel(C,1), Ce = gel(C,2), z = fracgammaeval(F,s,prec);
409 21280 : long i, lR = lg(Rw), lC = lg(Cw);
410 21280 : *ext = 0;
411 32214 : for (i = 1; i < lR; i++)
412 10934 : z = gmul(z, gpowgs(gamma_R(gadd(s,gel(Rw, i)), ext, prec), Re[i]));
413 32746 : for (i = 1; i < lC; i++)
414 11466 : z = gmul(z, gpowgs(gamma_C(gadd(s,gel(Cw, i)), ext, prec), Ce[i]));
415 21280 : return gc_upto(av, z);
416 : }
417 :
418 : static int
419 5446 : gammaordinary(GEN Vga, GEN s)
420 : {
421 5446 : long i, d = lg(Vga)-1;
422 14399 : for (i = 1; i <= d; i++)
423 : {
424 9142 : GEN z = gadd(s, gel(Vga,i));
425 : long e;
426 9142 : if (gexpo(imag_i(z)) < -10)
427 : {
428 9002 : z = real_i(z);
429 9002 : if (gsigne(z) <= 0) { (void)grndtoi(z, &e); if (e < -10) return 0; }
430 : }
431 : }
432 5257 : return 1;
433 : }
434 :
435 : /* Exponent A of t in asymptotic expansion; K(t) ~ C t^A exp(-pi d t^(2/d)).
436 : * suma = vecsum(Vga)*/
437 : static double
438 93165 : gammavec_expo(long d, double suma) { return (1 - d + suma) / d; }
439 :
440 : /*******************************************************************/
441 : /* First part: computations only involving Theta(t) */
442 : /*******************************************************************/
443 :
444 : static void
445 142072 : get_cone(GEN t, double *r, double *a)
446 : {
447 142072 : const long prec = LOWDEFAULTPREC;
448 142072 : if (typ(t) == t_COMPLEX)
449 : {
450 22078 : t = gprec_w(t, prec);
451 22078 : *r = gtodouble(gabs(t, prec));
452 22078 : *a = fabs(gtodouble(garg(t, prec)));
453 : }
454 : else
455 : {
456 119994 : *r = fabs(gtodouble(t));
457 119994 : *a = 0.;
458 : }
459 142072 : if (!*r && !*a) pari_err_DOMAIN("lfunthetainit","t","=",gen_0,t);
460 142065 : }
461 : /* slightly larger cone than necessary, to avoid round-off problems */
462 : static void
463 87318 : get_cone_fuzz(GEN t, double *r, double *a)
464 87318 : { get_cone(t, r, a); *r -= 1e-10; if (*a) *a += 1e-10; }
465 :
466 : /* Initialization m-th Theta derivative. tdom is either
467 : * - [rho,alpha]: assume |t| >= rho and |arg(t)| <= alpha
468 : * - a positive real scalar: assume t real, t >= tdom;
469 : * - a complex number t: compute at t;
470 : * N is the conductor (either the true one from ldata or a guess from
471 : * lfunconductor) */
472 : long
473 65555 : lfunthetacost(GEN ldata, GEN tdom, long m, long bit, long *extrabit)
474 : {
475 65555 : pari_sp av = avma;
476 65555 : GEN Vga = ldata_get_gammavec(ldata);
477 65555 : long d = lg(Vga)-1;
478 65555 : double k1 = maxdd(ldata_get_k1_dbl(ldata), 0.);
479 65555 : double c = d/2., a, A, B, logC, al, rho, T;
480 65555 : double N = gtodouble(ldata_get_conductor(ldata));
481 :
482 65555 : if (extrabit) *extrabit = 0;
483 65555 : if (!N) pari_err_TYPE("lfunthetaneed [missing conductor]", ldata);
484 65555 : if (typ(tdom) == t_VEC && lg(tdom) == 3)
485 : {
486 7 : rho= gtodouble(gel(tdom,1));
487 7 : al = gtodouble(gel(tdom,2));
488 : }
489 : else
490 65548 : get_cone_fuzz(tdom, &rho, &al);
491 65548 : A = gammavec_expo(d, gtodouble(sumVga(Vga))); set_avma(av);
492 65548 : a = (A+k1+1) + (m-1)/c;
493 65548 : if (fabs(a) < 1e-10) a = 0.;
494 65548 : logC = c*M_LN2 - log(c)/2;
495 : /* +1: fudge factor */
496 65548 : B = M_LN2*bit+logC+m*log(2*M_PI) + 1 + (k1+1)*log(N)/2 - (k1+m+1)*log(rho);
497 65548 : if (al)
498 : { /* t = rho e^(i*al), T^(1/c) = Re(t^(1/c)) > 0, T = rho cos^c(al/c) */
499 11046 : double z = cos(al/c);
500 11046 : if (z <= 0)
501 7 : pari_err_DOMAIN("lfunthetaneed", "arg t", ">", dbltor(c*M_PI/2), tdom);
502 11039 : T = (d == 2 && typ(tdom) != t_VEC)? gtodouble(real_i(tdom)): rho*pow(z,c);
503 11039 : B -= log(z) * (c * (k1+A+1) + m);
504 : }
505 : else
506 54502 : T = rho;
507 65541 : if (B <= 0) return 0;
508 65541 : A = floor(0.9 + dblcoro526(a,c,B) / T * sqrt(N));
509 65541 : if (dblexpo(A) >= BITS_IN_LONG-1) pari_err_OVERFLOW("lfunthetacost");
510 65534 : if (extrabit && A && a * log2(A) > 16) *extrabit = a * log2(A);
511 65534 : return (long)A;
512 : }
513 : long
514 21 : lfunthetacost0(GEN L, GEN tdom, long m, long bitprec)
515 : {
516 : long n;
517 21 : if (is_linit(L) && linit_get_type(L)==t_LDESC_THETA)
518 7 : {
519 7 : GEN tech = linit_get_tech(L);
520 7 : n = lg(theta_get_an(tech))-1;
521 : }
522 : else
523 : {
524 14 : pari_sp av = avma;
525 14 : GEN ldata = lfunmisc_to_ldata_shallow(L);
526 14 : n = lfunthetacost(ldata, tdom? tdom: gen_1, m, bitprec, NULL);
527 7 : set_avma(av);
528 : }
529 14 : return n;
530 : }
531 :
532 : static long
533 6923 : fracgammadegree(GEN FVga)
534 6923 : { GEN F = gel(FVga,1); return (typ(F)==t_RFRAC)? degpol(gel(F,2)): 0; }
535 :
536 : /* Poles of a L-function can be represented in the following ways:
537 : * 1) Nothing (ldata has only 6 components, ldata_get_residue = NULL).
538 : * 2) a complex number (single pole at s = k with given residue, unknown if 0).
539 : * 3) A vector (possibly empty) of 2-component vectors [a, ra], where a is the
540 : * pole, ra a t_SER: its Taylor expansion at a. A t_VEC encodes the polar
541 : * part of L, a t_COL, the polar part of Lambda */
542 :
543 : /* 'a' a complex number (pole), 'r' the polar part of L at 'a';
544 : * return 'R' the polar part of Lambda at 'a' */
545 : static GEN
546 5208 : rtoR(GEN a, GEN r, GEN FVga, GEN N, long prec)
547 : {
548 5208 : long v = lg(r)-2, d = fracgammadegree(FVga), ext;
549 5208 : GEN Na, as = deg1ser_shallow(gen_1, a, varn(r), v);
550 5208 : Na = gpow(N, gdivgu(as, 2), prec);
551 : /* make up for a possible loss of accuracy */
552 5208 : if (d) as = deg1ser_shallow(gen_1, a, varn(r), v + d);
553 5208 : return gmul(gmul(r, Na), gammafactproduct(FVga, as, &ext, prec));
554 : }
555 :
556 : /* assume r in normalized form: t_VEC of pairs [be,re] */
557 : GEN
558 4809 : lfunrtopoles(GEN r)
559 : {
560 4809 : long j, l = lg(r);
561 4809 : GEN v = cgetg(l, t_VEC);
562 10017 : for (j = 1; j < l; j++) gel(v,j) = gmael(r, j, 1);
563 4809 : gen_sort_inplace(v, (void*)&cmp_universal, cmp_nodata, NULL);
564 4809 : return v;
565 : }
566 :
567 : /* r / x + O(1) */
568 : static GEN
569 5236 : simple_pole(GEN r)
570 5236 : { return isintzero(r)? gen_0: serpole(r); }
571 : static GEN
572 6230 : normalize_simple_pole(GEN r, GEN k)
573 : {
574 6230 : long tx = typ(r);
575 6230 : if (is_vec_t(tx)) return r;
576 5236 : if (!is_scalar_t(tx)) pari_err_TYPE("lfunrootres [poles]", r);
577 5236 : return mkvec(mkvec2(k, simple_pole(r)));
578 : }
579 : /* check and normalize the description of a polar part as a t_VEC (r for L)
580 : * or t_COL (R for Lambda) of pairs [a, ra] */
581 : static GEN
582 5712 : normalizepoles(GEN r, GEN k)
583 : {
584 : long iv, j, l;
585 : GEN v;
586 5712 : if (!is_vec_t(typ(r))) return normalize_simple_pole(r, k);
587 2569 : v = cgetg_copy(r, &l);
588 6475 : for (j = iv = 1; j < l; j++)
589 : {
590 3906 : GEN rj = gel(r,j), a = gel(rj,1), ra = gel(rj,2);
591 3906 : if (!is_scalar_t(typ(a)) || typ(ra) != t_SER)
592 0 : pari_err_TYPE("lfunrootres [poles]",r);
593 3906 : gel(v,iv++) = rj;
594 : }
595 2569 : setlg(v, iv); return v;
596 : }
597 : static int
598 9212 : residues_known(GEN r)
599 : {
600 9212 : long i, l = lg(r);
601 9212 : if (isintzero(r)) return 0;
602 8883 : if (!is_vec_t(typ(r))) return 1;
603 12348 : for (i = 1; i < l; i++)
604 : {
605 7518 : GEN ri = gel(r,i);
606 7518 : if (!is_vec_t(typ(ri)) || lg(ri)!=3)
607 0 : pari_err_TYPE("lfunrootres [poles]",r);
608 7518 : if (isintzero(gel(ri, 2))) return 0;
609 : }
610 4830 : return 1;
611 : }
612 :
613 : /* Compute R's from r's (r = Taylor devts of L(s), R of Lambda(s)).
614 : * 'r/eno' passed to override the one from ldata */
615 : static GEN
616 24528 : lfunrtoR_i(GEN ldata, GEN r, GEN eno, long prec)
617 : {
618 24528 : GEN Vga = ldata_get_gammavec(ldata), N = ldata_get_conductor(ldata);
619 : GEN R, vr, FVga;
620 24528 : pari_sp av = avma;
621 : long lr, j, jR;
622 24528 : GEN k = ldata_get_k(ldata);
623 :
624 24528 : if (!r || isintzero(eno) || !residues_known(r))
625 18816 : return gen_0;
626 5712 : r = normalizepoles(r, k);
627 5712 : if (typ(r) == t_COL) return gc_GEN(av, r);
628 4809 : if (typ(ldata_get_dual(ldata)) != t_INT)
629 0 : pari_err(e_MISC,"please give the Taylor expansion of Lambda");
630 4809 : vr = lfunrtopoles(r); lr = lg(vr);
631 4809 : FVga = gammafactor(Vga);
632 4809 : R = cgetg(2*lr, t_COL);
633 10017 : for (j = jR = 1; j < lr; j++)
634 : {
635 5208 : GEN rj = gel(r,j), a = gel(rj,1), ra = gel(rj,2);
636 5208 : GEN b, Ra = rtoR(a, ra, FVga, N, prec);
637 5208 : long o = -valser(Ra); /* Lambda pole order */
638 5208 : if (o <= 0) continue;
639 5208 : b = gsub(k, conj_i(a));
640 5208 : if (lg(Ra)-2 < o)
641 0 : pari_err(e_MISC,
642 : "please give more terms in L function's Taylor expansion at %Ps", a);
643 5208 : setlg(Ra, o + 2); /* truncate useless terms */
644 5208 : gel(R,jR++) = mkvec2(a, Ra);
645 5208 : if (!tablesearch(vr, b, (int (*)(GEN,GEN))&cmp_universal))
646 : {
647 4991 : GEN mX = gneg(pol_x(varn(Ra)));
648 4991 : GEN Rb = gmul(eno, gsubst(conj_i(Ra), varn(Ra), mX));
649 4991 : gel(R,jR++) = mkvec2(b, Rb);
650 : }
651 : }
652 4809 : setlg(R, jR); return gc_GEN(av, R);
653 : }
654 : static GEN
655 24052 : lfunrtoR_eno(GEN ldata, GEN eno, long prec)
656 24052 : { return lfunrtoR_i(ldata, ldata_get_residue(ldata), eno, prec); }
657 : static GEN
658 21777 : lfunrtoR(GEN ldata, long prec)
659 21777 : { return lfunrtoR_eno(ldata, ldata_get_rootno(ldata), prec); }
660 :
661 : static long
662 21777 : prec_fix(long prec)
663 : {
664 : #ifndef LONG_IS_64BIT
665 : /* make sure that default accuracy is the same on 32/64bit */
666 3111 : if (odd(prec)) prec += EXTRAPREC64;
667 : #endif
668 21777 : return prec;
669 : }
670 :
671 : /* thetainit using {an: n <= L}; if (m = 0 && easytheta), an2 is an * n^al */
672 : static GEN
673 21777 : lfunthetainit0(GEN ldata, GEN tdom, GEN an2, long m,
674 : long bitprec, long extrabit)
675 : {
676 21777 : long prec = nbits2prec(bitprec);
677 21777 : GEN tech, N = ldata_get_conductor(ldata);
678 21777 : GEN K = gammamellininvinit(ldata, m, bitprec + extrabit);
679 21777 : GEN R = lfunrtoR(ldata, prec);
680 21777 : if (!tdom) tdom = gen_1;
681 21777 : if (typ(tdom) != t_VEC)
682 : {
683 : double r, a;
684 21770 : get_cone_fuzz(tdom, &r, &a);
685 21770 : tdom = mkvec2(dbltor(r), a? dbltor(a): gen_0);
686 : }
687 21777 : prec += maxss(EXTRAPREC64, nbits2extraprec(extrabit));
688 21777 : tech = mkvecn(7, an2,K,R, stoi(bitprec), stoi(m), tdom,
689 : gsqrt(ginv(N), prec_fix(prec)));
690 21777 : return mkvec3(mkvecsmall(t_LDESC_THETA), ldata, tech);
691 : }
692 :
693 : /* tdom: 1) positive real number r, t real, t >= r; or
694 : * 2) [r,a], describing the cone |t| >= r, |arg(t)| <= a */
695 : static GEN
696 10535 : lfunthetainit_i(GEN data, GEN tdom, long m, long bit)
697 : {
698 10535 : GEN ldata = lfunmisc_to_ldata_shallow(data);
699 10535 : long extrabit, b = 32, L = lfunthetacost(ldata, tdom, m, bit, &extrabit);
700 10521 : long prec = nbits2prec(bit + extrabit);
701 10521 : GEN ldatan = ldata_newprec(ldata, prec);
702 10521 : GEN an = ldata_get_an(ldatan), Vga = ldata_get_gammavec(ldatan);
703 10521 : an = ldata_vecan(an, L, prec);
704 10521 : if (m == 0 && Vgaeasytheta(Vga)) an = antwist(an, Vga, prec);
705 10521 : if (typ(an) != t_VECSMALL) b = maxss(b, gexpo(an));
706 10521 : return lfunthetainit0(ldatan, tdom, an, m, bit, b);
707 : }
708 :
709 : GEN
710 357 : lfunthetainit(GEN ldata, GEN tdom, long m, long bitprec)
711 : {
712 357 : pari_sp av = avma;
713 357 : GEN S = lfunthetainit_i(ldata, tdom? tdom: gen_1, m, bitprec);
714 357 : return gc_GEN(av, S);
715 : }
716 :
717 : GEN
718 2478 : lfunan(GEN ldata, long L, long prec)
719 : {
720 2478 : pari_sp av = avma;
721 : GEN an ;
722 2478 : ldata = ldata_newprec(lfunmisc_to_ldata_shallow(ldata), prec);
723 2478 : an = gc_GEN(av, ldata_vecan(ldata_get_an(ldata), L, prec));
724 2422 : if (typ(an) != t_VEC) an = vecsmall_to_vec_inplace(an);
725 2422 : return an;
726 : }
727 :
728 : static GEN
729 15890 : mulrealvec(GEN x, GEN y)
730 : {
731 15890 : if (is_vec_t(typ(x)) && is_vec_t(typ(y)))
732 84 : pari_APPLY_same(mulreal(gel(x,i),gel(y,i)))
733 : else
734 15862 : return mulreal(x,y);
735 : }
736 : static GEN
737 32114 : gmulvec(GEN x, GEN y)
738 : {
739 32114 : if (is_vec_t(typ(x)) && is_vec_t(typ(y)))
740 2702 : pari_APPLY_same(gmul(gel(x,i),gel(y,i)))
741 : else
742 31449 : return gmul(x,y);
743 : }
744 : static GEN
745 10129 : gdivvec(GEN x, GEN y)
746 : {
747 10129 : if (is_vec_t(typ(x)) && is_vec_t(typ(y)))
748 2247 : pari_APPLY_same(gdiv(gel(x,i),gel(y,i)))
749 : else
750 9555 : return gdiv(x,y);
751 : }
752 :
753 : static GEN
754 3647 : gsubvec(GEN x, GEN y)
755 : {
756 3647 : if (is_vec_t(typ(x)) && !is_vec_t(typ(y)))
757 0 : pari_APPLY_same(gsub(gel(x,i),y))
758 : else
759 3647 : return gsub(x,y);
760 : }
761 :
762 : /* return [1^(2/d), 2^(2/d),...,lim^(2/d)] */
763 : static GEN
764 9002 : mkvroots(long d, long lim, long prec)
765 : {
766 9002 : if (d <= 4)
767 : {
768 8652 : GEN v = cgetg(lim+1,t_VEC);
769 : long n;
770 8652 : switch(d)
771 : {
772 2870 : case 1:
773 53356 : for (n=1; n <= lim; n++) gel(v,n) = sqru(n);
774 2870 : return v;
775 1617 : case 2:
776 239197 : for (n=1; n <= lim; n++) gel(v,n) = utoipos(n);
777 1617 : return v;
778 2697 : case 4:
779 6969654 : for (n=1; n <= lim; n++) gel(v,n) = sqrtr(utor(n, prec));
780 2697 : return v;
781 : }
782 : }
783 1818 : return vecpowug(lim, gdivgu(gen_2,d), prec);
784 : }
785 :
786 : GEN
787 62965 : lfunthetacheckinit(GEN data, GEN t, long m, long bitprec)
788 : {
789 62965 : if (is_linit(data) && linit_get_type(data)==t_LDESC_THETA)
790 : {
791 54754 : GEN tdom, thetainit = linit_get_tech(data);
792 54754 : long bitprecnew = theta_get_bitprec(thetainit);
793 54754 : long m0 = theta_get_m(thetainit);
794 : double r, al, rt, alt;
795 54754 : if (m0 != m)
796 0 : pari_err_DOMAIN("lfuntheta","derivative order","!=", stoi(m),stoi(m0));
797 54754 : if (bitprec > bitprecnew) goto INIT;
798 54754 : get_cone(t, &rt, &alt);
799 54754 : tdom = theta_get_tdom(thetainit);
800 54754 : r = gtodouble(gel(tdom,1));
801 54754 : al= gtodouble(gel(tdom,2)); if (rt >= r && alt <= al) return data;
802 : }
803 8211 : INIT:
804 10087 : return lfunthetainit_i(data, t, m, bitprec);
805 : }
806 :
807 : static GEN
808 15159549 : get_an(GEN an, long n)
809 : {
810 15159549 : if (typ(an) == t_VECSMALL) { long a = an[n]; if (a) return stoi(a); }
811 15159549 : else { GEN a = gel(an,n); if (a && !gequal0(a)) return a; }
812 12822350 : return NULL;
813 : }
814 : /* x * an[n] */
815 : static GEN
816 14693361 : mul_an(GEN an, long n, GEN x)
817 : {
818 14693361 : if (typ(an) == t_VECSMALL) { long a = an[n]; if (a) return gmulsg(a,x); }
819 9431229 : else { GEN a = gel(an,n); if (a && !gequal0(a)) return gmul(a,x); }
820 2611979 : return NULL;
821 : }
822 : /* 2*t^a * x **/
823 : static GEN
824 334733 : mulT(GEN t, GEN a, GEN x, long prec)
825 : {
826 334733 : if (gequal0(a)) return gmul2n(x,1);
827 32725 : return gmul(x, gmul2n(gequal1(a)? t: gpow(t,a,prec), 1));
828 : }
829 :
830 : static GEN
831 35699737 : vecan_cmul(void *E, GEN P, long a, GEN x)
832 : {
833 : (void)E;
834 35699737 : if (typ(P) == t_VECSMALL)
835 24762599 : return (a==0 || !P[a])? NULL: gmulsg(P[a], x);
836 : else
837 10937138 : return (a==0 || !gel(P,a))? NULL: gmul(gel(P,a), x);
838 : }
839 : /* d=2, 2 sum_{n <= N} a(n) (n t)^al q^n, q = exp(-2pi t),
840 : * an2[n] = a(n) * n^al */
841 : static GEN
842 295095 : theta2_i(GEN an2, long N, GEN t, GEN al, long prec)
843 : {
844 295095 : GEN S, q, pi2 = Pi2n(1,prec);
845 295095 : const struct bb_algebra *alg = get_Rg_algebra();
846 295095 : setsigne(pi2,-1); q = gexp(gmul(pi2, t), prec);
847 : /* Brent-Kung in case the a_n are small integers */
848 295095 : S = gen_bkeval(an2, N, q, 1, NULL, alg, vecan_cmul);
849 295095 : return mulT(t, al, S, prec);
850 : }
851 : static GEN
852 286923 : theta2(GEN an2, long N, GEN t, GEN al, long prec)
853 : {
854 286923 : pari_sp av = avma;
855 286923 : return gc_upto(av, theta2_i(an2, N, t, al, prec));
856 : }
857 :
858 : /* d=1, 2 sum_{n <= N} a_n (n t)^al q^(n^2), q = exp(-pi t^2),
859 : * an2[n] is a_n n^al */
860 : static GEN
861 39638 : theta1(GEN an2, long N, GEN t, GEN al, long prec)
862 : {
863 39638 : GEN q = gexp(gmul(negr(mppi(prec)), gsqr(t)), prec);
864 39638 : GEN vexp = gsqrpowers(q, N), S = gen_0;
865 39638 : pari_sp av = avma;
866 : long n;
867 6621115 : for (n = 1; n <= N; n++)
868 : {
869 6581477 : GEN c = mul_an(an2, n, gel(vexp,n));
870 6581477 : if (c)
871 : {
872 5505067 : S = gadd(S, c);
873 5505067 : if (gc_needed(av, 3)) S = gc_upto(av, S);
874 : }
875 : }
876 39638 : return mulT(t, al, S, prec);
877 : }
878 :
879 : /* If m > 0, compute m-th derivative of theta(t) = theta0(t/sqrt(N))
880 : * with absolute error 2^-bitprec; theta(t)=\sum_{n\ge1}a(n)K(nt/N^(1/2)) */
881 : GEN
882 52486 : lfuntheta(GEN data, GEN t, long m, long bitprec)
883 : {
884 52486 : pari_sp ltop = avma;
885 : long limt, d;
886 : GEN isqN, vecan, Vga, ldata, theta, thetainit, S;
887 : long n, prec;
888 :
889 52486 : theta = lfunthetacheckinit(data, t, m, bitprec);
890 52479 : ldata = linit_get_ldata(theta);
891 52479 : thetainit = linit_get_tech(theta);
892 52479 : vecan = theta_get_an(thetainit);
893 52479 : isqN = theta_get_isqrtN(thetainit);
894 52479 : prec = maxss(realprec(isqN), nbits2prec(bitprec));
895 52479 : t = gprec_w(t, prec);
896 52479 : limt = lg(vecan)-1;
897 52479 : if (theta == data)
898 48265 : limt = minss(limt, lfunthetacost(ldata, t, m, bitprec, NULL));
899 52479 : if (!limt)
900 : {
901 14 : set_avma(ltop); S = real_0_bit(-bitprec);
902 14 : if (!is_real_t(typ(t)) || !ldata_isreal(ldata))
903 7 : S = gc_GEN(ltop, mkcomplex(S,S));
904 14 : return S;
905 : }
906 52465 : t = gmul(t, isqN);
907 52465 : Vga = ldata_get_gammavec(ldata);
908 52465 : d = lg(Vga)-1;
909 52465 : if (m == 0 && Vgaeasytheta(Vga))
910 : {
911 47810 : if (theta_get_m(thetainit) > 0) vecan = antwist(vecan, Vga, prec);
912 47810 : if (d == 1) S = theta1(vecan, limt, t, gel(Vga,1), prec);
913 8172 : else S = theta2_i(vecan, limt, t, vecmin(Vga), prec);
914 : }
915 : else
916 : {
917 4655 : GEN K = theta_get_K(thetainit);
918 4655 : GEN vroots = mkvroots(d, limt, prec);
919 : pari_sp av;
920 4655 : t = gpow(t, gdivgu(gen_2,d), prec);
921 4655 : S = gen_0; av = avma;
922 15164204 : for (n = 1; n <= limt; ++n)
923 : {
924 15159549 : GEN nt, an = get_an(vecan, n);
925 15159549 : if (!an) continue;
926 2337199 : nt = gmul(gel(vroots,n), t);
927 2337199 : if (m) an = gmul(an, powuu(n, m));
928 2337199 : S = gadd(S, gmul(an, gammamellininvrt(K, nt, bitprec)));
929 2337199 : if ((n & 0x1ff) == 0) S = gc_upto(av, S);
930 : }
931 4655 : if (m) S = gmul(S, gpowgs(isqN, m));
932 : }
933 52465 : return gc_upto(ltop, S);
934 : }
935 :
936 : /*******************************************************************/
937 : /* Second part: Computation of L-Functions. */
938 : /*******************************************************************/
939 :
940 : struct lfunp {
941 : long precmax, Dmax, D, M, m0, nmax, d, vgaell;
942 : double k1, dc, dw, dh, MAXs, sub;
943 : GEN L, an, bn;
944 : };
945 :
946 : static void
947 27617 : lfunp_set(GEN ldata, long der, long bitprec, struct lfunp *S)
948 : {
949 27617 : const long derprec = (der > 1)? dbllog2(mpfact(der)): 0; /* log2(der!) */
950 : GEN Vga, N, L, k;
951 : long k1, d, m, M, flag, nmax;
952 : double a, A, E, hd, Ep, d2, suma, maxs, mins, sub, B0,B1;
953 : double logN2, logC, Lestimate, Mestimate;
954 :
955 27617 : Vga = ldata_get_gammavec(ldata);
956 27617 : S->d = d = lg(Vga)-1; d2 = d/2.;
957 :
958 27617 : suma = gtodouble(sumVga(Vga));
959 27617 : k = ldata_get_k(ldata);
960 27617 : N = ldata_get_conductor(ldata);
961 27617 : logN2 = log(gtodouble(N)) / 2;
962 27617 : maxs = S->dc + S->dw;
963 27617 : mins = S->dc - S->dw;
964 27617 : S->MAXs = maxdd(maxs, gtodouble(k)-mins);
965 :
966 : /* we compute Lambda^(der)(s) / der!; need to compensate for L^(der)(s)
967 : * ln |gamma(s)| ~ -(pi/4) \sum_i |Im(s + a_i)|; max with 1: fudge factor */
968 27617 : a = (M_PI/(4*M_LN2))*(d*S->dh + sumVgaimpos(Vga));
969 27617 : S->D = (long)ceil(bitprec + derprec + maxdd(a, 1));
970 27617 : E = M_LN2*S->D; /* D:= required absolute bitprec */
971 :
972 27617 : Ep = E + maxdd(M_PI * S->dh * d2, (d*S->MAXs + suma - 1) * log(E));
973 27617 : hd = d2*M_PI*M_PI / Ep;
974 27617 : S->m0 = (long)ceil(M_LN2/hd);
975 27617 : hd = M_LN2/S->m0;
976 :
977 27617 : logC = d2*M_LN2 - log(d2)/2;
978 27617 : k1 = maxdd(ldata_get_k1_dbl(ldata), 0.);
979 27617 : S->k1 = k1; /* assume |a_n| << n^k1 with small implied constant */
980 27617 : A = gammavec_expo(d, suma);
981 :
982 27617 : sub = 0.;
983 27617 : if (mins > 1)
984 : {
985 5446 : GEN sig = dbltor(mins);
986 5446 : sub += logN2*mins;
987 5446 : if (gammaordinary(Vga, sig))
988 : {
989 : long ext;
990 5257 : GEN gas = gammafactproduct(gammafactor(Vga), sig, &ext, LOWDEFAULTPREC);
991 5257 : if (typ(gas) != t_SER)
992 : {
993 5257 : double dg = dbllog2(gas);
994 5257 : if (dg > 0) sub += dg * M_LN2;
995 : }
996 : }
997 : }
998 27617 : S->sub = sub;
999 27617 : M = 1000;
1000 27617 : L = cgetg(M+2, t_VECSMALL);
1001 27617 : a = S->k1 + A;
1002 :
1003 27617 : B0 = 5 + E - S->sub + logC + S->k1*logN2; /* 5 extra bits */
1004 27617 : B1 = hd * (S->MAXs - S->k1);
1005 27617 : Lestimate = dblcoro526(a + S->MAXs - 2./d, d/2.,
1006 27617 : E - S->sub + logC - log(2*M_PI*hd) + S->MAXs*logN2);
1007 27617 : Mestimate = ((Lestimate > 0? log(Lestimate): 0) + logN2) / hd;
1008 27617 : nmax = 0;
1009 27617 : flag = 0;
1010 27617 : for (m = 0;; m++)
1011 2439332 : {
1012 2466949 : double x, H = logN2 - m*hd, B = B0 + m*B1;
1013 : long n;
1014 2466949 : x = dblcoro526(a, d/2., B);
1015 2466949 : n = floor(x*exp(H));
1016 2466949 : if (n > nmax) nmax = n;
1017 2466949 : if (m > M) { M *= 2; L = vecsmall_lengthen(L,M+2); }
1018 2466949 : L[m+1] = n;
1019 2466949 : if (n == 0) { if (++flag > 2 && m > Mestimate) break; } else flag = 0;
1020 : }
1021 28457 : m -= 2; while (m > 0 && !L[m]) m--;
1022 27617 : if (m == 0) { nmax = 1; L[1] = 1; m = 1; } /* can happen for tiny bitprec */
1023 27617 : setlg(L, m+1); S->M = m-1;
1024 27617 : S->L = L;
1025 27617 : S->nmax = nmax;
1026 :
1027 27617 : S->Dmax = S->D + (long)ceil((S->M * hd * S->MAXs - S->sub) / M_LN2);
1028 27617 : if (S->Dmax < S->D) S->Dmax = S->D;
1029 27617 : S->precmax = nbits2prec(S->Dmax);
1030 27617 : if (DEBUGLEVEL > 1)
1031 0 : err_printf("Dmax=%ld, D=%ld, M = %ld, nmax = %ld, m0 = %ld\n",
1032 : S->Dmax,S->D,S->M,S->nmax, S->m0);
1033 27617 : }
1034 :
1035 : static GEN
1036 11410 : lfuninit_pol(GEN v, GEN poqk, long prec)
1037 : {
1038 11410 : long m, M = lg(v) - 2;
1039 11410 : GEN pol = cgetg(M+3, t_POL);
1040 11410 : pol[1] = evalsigne(1) | evalvarn(0);
1041 11410 : gel(pol, 2) = gprec_w(gmul2n(gel(v,1), -1), prec);
1042 11410 : if (poqk)
1043 544988 : for (m = 2; m <= M+1; m++)
1044 533634 : gel(pol, m+1) = gprec_w(gmul(gel(poqk,m), gel(v,m)), prec);
1045 : else
1046 2324 : for (m = 2; m <= M+1; m++)
1047 2268 : gel(pol, m+1) = gprec_w(gel(v,m), prec);
1048 11410 : return RgX_renormalize_lg(pol, M+3);
1049 : }
1050 :
1051 : static void
1052 81995 : worker_init(long q, GEN *an, GEN *bn, GEN *AB, GEN *A, GEN *B)
1053 : {
1054 81995 : if (typ(*bn) == t_INT) *bn = NULL;
1055 81995 : if (*bn)
1056 : {
1057 712 : *AB = cgetg(3, t_VEC);
1058 712 : gel(*AB,1) = *A = cgetg(q+1, t_VEC);
1059 712 : gel(*AB,2) = *B = cgetg(q+1, t_VEC);
1060 712 : if (typ(an) == t_VEC) *an = RgV_kill0(*an);
1061 712 : if (typ(bn) == t_VEC) *bn = RgV_kill0(*bn);
1062 : }
1063 : else
1064 : {
1065 81283 : *B = NULL;
1066 81283 : *AB = *A = cgetg(q+1, t_VEC);
1067 81283 : if (typ(*an) == t_VEC) *an = RgV_kill0(*an);
1068 : }
1069 81995 : }
1070 : GEN
1071 22396 : lfuninit_theta2_worker(long r, GEN L, GEN qk, GEN a, GEN di, GEN an, GEN bn)
1072 : {
1073 22396 : long q, m, prec = di[1], M = di[2], m0 = di[3], L0 = lg(an)-1;
1074 : GEN AB, A, B;
1075 22396 : worker_init((M - r) / m0 + 1, &an, &bn, &AB, &A, &B);
1076 302354 : for (q = 0, m = r; m <= M; m += m0, q++)
1077 : {
1078 279958 : GEN t = gel(qk, m+1);
1079 279958 : long N = minss(L[m+1],L0);
1080 279958 : gel(A, q+1) = theta2(an, N, t, a, prec); /* theta(exp(mh)) */
1081 279958 : if (bn) gel(B, q+1) = theta2(bn, N, t, a, prec);
1082 : }
1083 22396 : return AB;
1084 : }
1085 :
1086 : /* theta(exp(mh)) ~ sum_{n <= N} a(n) k[m,n] */
1087 : static GEN
1088 260389 : an_msum(GEN an, long N, GEN vKm)
1089 : {
1090 260389 : pari_sp av = avma;
1091 260389 : GEN s = gen_0;
1092 : long n;
1093 15140909 : for (n = 1; n <= N; n++)
1094 14880520 : if (gel(vKm,n))
1095 : {
1096 8111884 : GEN c = mul_an(an, n, gel(vKm,n));
1097 8111884 : if (c) s = gadd(s, c);
1098 : }
1099 260389 : return gc_upto(av, s);
1100 : }
1101 :
1102 : GEN
1103 59599 : lfuninit_worker(long r, GEN K, GEN L, GEN peh2d, GEN vroots, GEN dr, GEN di,
1104 : GEN an, GEN bn)
1105 : {
1106 59599 : pari_sp av0 = avma;
1107 59599 : long m, n, q, L0 = lg(an)-1;
1108 59599 : double sig0 = rtodbl(gel(dr,1)), sub2 = rtodbl(gel(dr,2));
1109 59599 : double k1 = rtodbl(gel(dr,3)), MAXs = rtodbl(gel(dr,4));
1110 59599 : long D = di[1], M = di[2], m0 = di[3];
1111 59599 : double M0 = sig0? sub2 / sig0: 1./0.;
1112 59599 : GEN AB, A, B, vK = cgetg(M/m0 + 2, t_VEC);
1113 :
1114 318973 : for (q = 0, m = r; m <= M; m += m0, q++)
1115 259374 : gel(vK, q+1) = const_vec(L[m+1], NULL);
1116 59599 : worker_init(q, &an, &bn, &AB, &A, &B);
1117 318973 : for (m -= m0, q--; m >= 0; m -= m0, q--)
1118 : {
1119 259374 : double c1 = D + ((m > M0)? m * sig0 - sub2 : 0);
1120 259374 : GEN vKm = gel(vK,q+1); /* conceptually K(m,n) */
1121 15135799 : for (n = 1; n <= L[m+1]; n++)
1122 : {
1123 : GEN t2d, kmn;
1124 14876425 : long nn, mm, qq, p = 0;
1125 : double c, c2;
1126 : pari_sp av;
1127 :
1128 14876425 : if (gel(vKm, n)) continue; /* done already */
1129 10806984 : c = c1 + k1 * log2(n);
1130 : /* n *= 2; m -= m0 => c += c2, provided m >= M0. Else c += k1 */
1131 10806984 : c2 = k1 - MAXs;
1132 : /* p = largest (absolute) accuracy to which we need K(m,n) */
1133 17682995 : for (mm=m,nn=n; mm >= M0;)
1134 : {
1135 13976155 : if (nn <= L[mm+1] && (gel(an, nn) || (bn && gel(bn, nn))))
1136 4731363 : if (c > 0) p = maxuu(p, (ulong)c);
1137 13976155 : nn <<= 1;
1138 13976155 : mm -= m0; if (mm >= M0) c += c2; else { c += k1; break; }
1139 : }
1140 : /* mm < M0 || nn > L[mm+1] */
1141 18216344 : for ( ; mm >= 0; nn<<=1,mm-=m0,c+=k1)
1142 7409360 : if (nn <= L[mm+1] && (gel(an, nn) || (bn && gel(bn, nn))))
1143 1841221 : if (c > 0) p = maxuu(p, (ulong)c);
1144 10806984 : if (!p) continue; /* a_{n 2^v} = 0 for all v in range */
1145 4038348 : av = avma;
1146 4038348 : t2d = mpmul(gel(vroots,n), gel(peh2d,m+1));/*(n exp(mh)/sqrt(N))^(2/d)*/
1147 4038348 : kmn = gc_upto(av, gammamellininvrt(K, t2d, p));
1148 12321809 : for (qq=q,mm=m,nn=n; mm >= 0; nn<<=1,mm-=m0,qq--)
1149 8283461 : if (nn <= L[mm+1]) gmael(vK, qq+1, nn) = kmn;
1150 : }
1151 : }
1152 318973 : for (q = 0, m = r; m <= M; m += m0, q++)
1153 : {
1154 259374 : long N = minss(L0, L[m+1]);
1155 259374 : gel(A, q+1) = an_msum(an, N, gel(vK,q+1));
1156 259374 : if (bn) gel(B, q+1) = an_msum(bn, N, gel(vK,q+1));
1157 : }
1158 59599 : return gc_upto(av0, AB);
1159 : }
1160 : /* return A = [\theta(exp(mh)), m=0..M], theta(t) = sum a(n) K(n/sqrt(N) t),
1161 : * h = log(2)/m0. If bn != NULL, return the pair [A, B] */
1162 : static GEN
1163 11256 : lfuninit_ab(GEN theta, GEN h, struct lfunp *S)
1164 : {
1165 11256 : const long M = S->M, prec = S->precmax;
1166 11256 : GEN tech = linit_get_tech(theta), isqN = theta_get_isqrtN(tech);
1167 11256 : GEN an = S->an, bn = S->bn, va, vb;
1168 : struct pari_mt pt;
1169 : GEN worker;
1170 : long m0, r, pending;
1171 :
1172 11256 : if (S->vgaell)
1173 : { /* d=2 and Vga = [a,a+1] */
1174 7126 : GEN a = vecmin(ldata_get_gammavec(linit_get_ldata(theta)));
1175 7126 : GEN qk = gpowers0(mpexp(h), M, isqN);
1176 7126 : m0 = minss(M+1, mt_nbthreads());
1177 7126 : worker = snm_closure(is_entry("_lfuninit_theta2_worker"),
1178 : mkvecn(6, S->L, qk, a, mkvecsmall3(prec, M, m0),
1179 : an, bn? bn: gen_0));
1180 : }
1181 : else
1182 : {
1183 : GEN vroots, peh2d, d2;
1184 4130 : double sig0 = S->MAXs / S->m0, sub2 = S->sub / M_LN2;
1185 : /* For all 0<= m <= M, and all n <= L[m+1] such that a_n!=0, we compute
1186 : * k[m,n] = K(n exp(mh)/sqrt(N))
1187 : * with ln(absolute error) <= E + max(mh sigma - sub, 0) + k1 * log(n).
1188 : * N.B. we use the 'rt' variant and pass (n exp(mh)/sqrt(N))^(2/d).
1189 : * Speedup: if n' = 2n and m' = m - m0 >= 0; then k[m,n] = k[m',n']. */
1190 4130 : vroots = mkvroots(S->d, S->nmax, prec); /* vroots[n] = n^(2/d) */
1191 4130 : d2 = gdivgu(gen_2, S->d);
1192 4130 : peh2d = gpowers0(gexp(gmul(d2,h), prec), M, gpow(isqN, d2, prec));
1193 4130 : m0 = S->m0; /* peh2d[m+1] = (exp(mh)/sqrt(N))^(2/d) */
1194 4130 : worker = snm_closure(is_entry("_lfuninit_worker"),
1195 : mkvecn(8, theta_get_K(tech), S->L, peh2d, vroots,
1196 : mkvec4(dbltor(sig0), dbltor(sub2),
1197 : dbltor(S->k1), dbltor(S->MAXs)),
1198 : mkvecsmall3(S->D, M, m0),
1199 : an, bn? bn: gen_0));
1200 : /* For each 0 <= m <= M, we will sum for n<=L[m+1] a(n) K(m,n)
1201 : * bit accuracy for K(m,n): D + k1*log2(n) + 1_{m > M0} (m*sig0 - sub2)
1202 : * We restrict m to arithmetic progressions r mod m0 to save memory and
1203 : * allow parallelization */
1204 : }
1205 11256 : va = cgetg(M+2, t_VEC);
1206 11256 : vb = bn? cgetg(M+2, t_VEC): NULL;
1207 11256 : mt_queue_start_lim(&pt, worker, m0);
1208 11256 : pending = 0;
1209 114861 : for (r = 0; r < m0 || pending; r++)
1210 : { /* m = q m0 + r */
1211 : GEN done, A, B;
1212 : long q, m, workid;
1213 103605 : mt_queue_submit(&pt, r, r < m0 ? mkvec(utoi(r)): NULL);
1214 103605 : done = mt_queue_get(&pt, &workid, &pending);
1215 103605 : if (!done) continue;
1216 81995 : if (bn) { A = gel(done,1); B = gel(done,2); } else { A = done; B = NULL; }
1217 621327 : for (q = 0, m = workid; m <= M; m += m0, q++)
1218 : {
1219 539332 : gel(va, m+1) = gel(A, q+1);
1220 539332 : if (bn) gel(vb, m+1) = gel(B, q+1);
1221 : }
1222 : }
1223 11256 : mt_queue_end(&pt);
1224 11256 : return bn? mkvec2(va, vb): va;
1225 : }
1226 :
1227 : static void
1228 143121 : parse_dom(double k, GEN dom, struct lfunp *S)
1229 : {
1230 143121 : long l = lg(dom);
1231 143121 : if (typ(dom)!=t_VEC) pari_err_TYPE("lfuninit [domain]", dom);
1232 143121 : if (l == 1)
1233 : {
1234 98 : S->dc = 0;
1235 98 : S->dw = -1;
1236 98 : S->dh = -1; return;
1237 : }
1238 143023 : if (l == 2)
1239 : {
1240 38110 : S->dc = k/2.;
1241 38110 : S->dw = 0.;
1242 38110 : S->dh = gtodouble(gel(dom,1));
1243 : }
1244 104913 : else if (l == 3)
1245 : {
1246 301 : S->dc = k/2.;
1247 301 : S->dw = gtodouble(gel(dom,1));
1248 301 : S->dh = gtodouble(gel(dom,2));
1249 : }
1250 104612 : else if (l == 4)
1251 : {
1252 104612 : S->dc = gtodouble(gel(dom,1));
1253 104612 : S->dw = gtodouble(gel(dom,2));
1254 104612 : S->dh = gtodouble(gel(dom,3));
1255 : }
1256 : else
1257 : {
1258 0 : pari_err_TYPE("lfuninit [domain]", dom);
1259 0 : S->dc = S->dw = S->dh = 0; /*-Wall*/
1260 : }
1261 143023 : if (S->dw < 0 || S->dh < 0) pari_err_TYPE("lfuninit [domain]", dom);
1262 : }
1263 :
1264 : /* do we have dom \subset dom0 ? dom = [center, width, height] */
1265 : int
1266 25190 : sdomain_isincl(double k, GEN dom, GEN dom0)
1267 : {
1268 : struct lfunp S0, S;
1269 25190 : parse_dom(k, dom, &S); if (S.dw < 0) return 1;
1270 25190 : parse_dom(k, dom0, &S0); if (S0.dw < 0) return 0;
1271 25190 : return S0.dc - S0.dw <= S.dc - S.dw
1272 25190 : && S0.dc + S0.dw >= S.dc + S.dw && S0.dh >= S.dh;
1273 : }
1274 :
1275 : static int
1276 25267 : checklfuninit(GEN linit, GEN DOM, long der, long bitprec)
1277 : {
1278 25267 : GEN ldata = linit_get_ldata(linit);
1279 25267 : GEN domain = lfun_get_domain(linit_get_tech(linit));
1280 25267 : GEN dom = domain_get_dom(domain);
1281 25267 : if (lg(dom) == 1) return 1;
1282 25190 : return domain_get_der(domain) >= der
1283 25190 : && domain_get_bitprec(domain) >= bitprec
1284 50380 : && sdomain_isincl(gtodouble(ldata_get_k(ldata)), DOM, dom);
1285 : }
1286 :
1287 : static GEN
1288 2394 : ginvsqrtvec(GEN x, long prec)
1289 : {
1290 2394 : if (is_vec_t(typ(x)))
1291 1813 : pari_APPLY_same(ginv(gsqrt(gel(x,i), prec)))
1292 1946 : else return ginv(gsqrt(x, prec));
1293 : }
1294 :
1295 : GEN
1296 12362 : lfuninit_make(long t, GEN ldata, GEN tech, GEN domain)
1297 : {
1298 12362 : GEN Vga = ldata_get_gammavec(ldata);
1299 12362 : long d = lg(Vga)-1;
1300 12362 : GEN w2 = gen_1, k2 = gmul2n(ldata_get_k(ldata), -1);
1301 12362 : GEN expot = gdivgu(gadd(gmulsg(d, gsubgs(k2, 1)), sumVga(Vga)), 4);
1302 12362 : if (typ(ldata_get_dual(ldata))==t_INT)
1303 : {
1304 12208 : GEN eno = ldata_get_rootno(ldata);
1305 12208 : long prec = nbits2prec( domain_get_bitprec(domain) );
1306 12208 : if (!isint1(eno)) w2 = ginvsqrtvec(eno, prec);
1307 : }
1308 12362 : tech = mkvec3(domain, tech, mkvec4(k2, w2, expot, gammafactor(Vga)));
1309 12362 : return mkvec3(mkvecsmall(t), ldata, tech);
1310 : }
1311 : static GEN
1312 224 : lfunnoinit(GEN ldata, long bitprec)
1313 : {
1314 224 : GEN tech, domain = mkvec2(cgetg(1,t_VEC), mkvecsmall2(0, bitprec));
1315 224 : GEN R = gen_0, r = ldata_get_residue(ldata), v = lfunrootres(ldata, bitprec);
1316 224 : ldata = shallowcopy(ldata);
1317 224 : gel(ldata,6) = gel(v,3);
1318 224 : if (r)
1319 : {
1320 196 : if (isintzero(r)) setlg(ldata,7); else gel(ldata,7) = r;
1321 196 : R = gel(v,2);
1322 : }
1323 224 : tech = mkvec3(domain, gen_0, R);
1324 224 : return lfuninit_make(t_LDESC_INIT, ldata, tech, domain);
1325 : }
1326 :
1327 : static void
1328 4130 : lfunparams2(struct lfunp *S)
1329 : {
1330 4130 : GEN L = S->L, an = S->an, bn = S->bn;
1331 : double pmax;
1332 4130 : long m, nan, nmax, neval, M = S->M;
1333 :
1334 4130 : S->vgaell = 0;
1335 : /* try to reduce parameters now we know the a_n (some may be 0) */
1336 4130 : if (typ(an) == t_VEC) an = RgV_kill0(an);
1337 4130 : nan = S->nmax; /* lg(an)-1 may be large than this */
1338 4130 : nmax = neval = 0;
1339 4130 : if (!bn)
1340 262468 : for (m = 0; m <= M; m++)
1341 : {
1342 258359 : long n = minss(nan, L[m+1]);
1343 369890 : while (n > 0 && !gel(an,n)) n--;
1344 258359 : if (n > nmax) nmax = n;
1345 258359 : neval += n;
1346 258359 : L[m+1] = n; /* reduce S->L[m+1] */
1347 : }
1348 : else
1349 : {
1350 21 : if (typ(bn) == t_VEC) bn = RgV_kill0(bn);
1351 1036 : for (m = 0; m <= M; m++)
1352 : {
1353 1015 : long n = minss(nan, L[m+1]);
1354 1057 : while (n > 0 && !gel(an,n) && !gel(bn,n)) n--;
1355 1015 : if (n > nmax) nmax = n;
1356 1015 : neval += n;
1357 1015 : L[m+1] = n; /* reduce S->L[m+1] */
1358 : }
1359 : }
1360 4130 : if (DEBUGLEVEL >= 1) err_printf("expected evaluations: %ld\n", neval);
1361 4130 : for (; M > 0; M--)
1362 4130 : if (L[M+1]) break;
1363 4130 : setlg(L, M+2);
1364 4130 : S->M = M;
1365 4130 : S->nmax = nmax;
1366 :
1367 : /* need K(n*exp(mh)/sqrt(N)) to absolute accuracy
1368 : * D + k1*log(n) + max(m * sig0 - sub2, 0) */
1369 4130 : pmax = S->D + S->k1 * log2(L[1]);
1370 4130 : if (S->MAXs)
1371 : {
1372 4130 : double sig0 = S->MAXs/S->m0, sub2 = S->sub / M_LN2;
1373 222896 : for (m = ceil(sub2 / sig0); m <= S->M; m++)
1374 : {
1375 218766 : double c = S->D + m*sig0 - sub2;
1376 218766 : if (S->k1 > 0) c += S->k1 * log2(L[m+1]);
1377 218766 : pmax = maxdd(pmax, c);
1378 : }
1379 : }
1380 4130 : S->Dmax = pmax;
1381 4130 : S->precmax = nbits2prec(pmax);
1382 4130 : }
1383 :
1384 : static GEN
1385 11270 : lfun_init_theta(GEN ldata, GEN eno, struct lfunp *S)
1386 : {
1387 11270 : GEN an2, dual, tdom = NULL, Vga = ldata_get_gammavec(ldata);
1388 11270 : long L, extrabit = 0, prec = S->precmax;
1389 11270 : if (eno)
1390 6650 : L = S->nmax;
1391 : else
1392 : {
1393 4620 : tdom = dbltor(sqrt(0.5));
1394 4620 : L = maxss(S->nmax, lfunthetacost(ldata, tdom, 0, S->D, &extrabit));
1395 4620 : prec += nbits2extraprec(extrabit);
1396 : }
1397 11270 : dual = ldata_get_dual(ldata);
1398 11270 : S->an = ldata_vecan(ldata_get_an(ldata), L, prec);
1399 11256 : S->bn = typ(dual)==t_INT? NULL: ldata_vecan(dual, S->nmax, prec);
1400 11256 : if (!vgaell(Vga)) lfunparams2(S);
1401 : else
1402 : {
1403 7126 : S->an = antwist(S->an, Vga, prec);
1404 7126 : if (S->bn) S->bn = antwist(S->bn, Vga, prec);
1405 7126 : S->vgaell = 1;
1406 : }
1407 11256 : an2 = lg(Vga)-1 == 1? antwist(S->an, Vga, prec): S->an;
1408 11256 : return lfunthetainit0(ldata, tdom, an2, 0, S->Dmax, extrabit);
1409 : }
1410 :
1411 : GEN
1412 16347 : lfuncost(GEN L, GEN dom, long der, long bit)
1413 : {
1414 16347 : pari_sp av = avma;
1415 16347 : GEN ldata = lfunmisc_to_ldata_shallow(L);
1416 16347 : GEN w, k = ldata_get_k(ldata);
1417 : struct lfunp S;
1418 :
1419 16347 : parse_dom(gtodouble(k), dom, &S); if (S.dw < 0) return mkvecsmall2(0, 0);
1420 16347 : lfunp_set(ldata, der, bit, &S);
1421 16347 : w = ldata_get_rootno(ldata);
1422 16347 : if (isintzero(w)) /* for lfunrootres */
1423 7 : S.nmax = maxss(S.nmax, lfunthetacost(ldata,dbltor(sqrt(0.5)),0,bit+1,NULL));
1424 16347 : set_avma(av); return mkvecsmall2(S.nmax, S.Dmax);
1425 : }
1426 : GEN
1427 49 : lfuncost0(GEN L, GEN dom, long der, long bitprec)
1428 : {
1429 49 : pari_sp av = avma;
1430 : GEN C;
1431 :
1432 49 : if (is_linit(L))
1433 : {
1434 28 : GEN tech = linit_get_tech(L);
1435 28 : GEN domain = lfun_get_domain(tech);
1436 28 : dom = domain_get_dom(domain);
1437 28 : der = domain_get_der(domain);
1438 28 : bitprec = domain_get_bitprec(domain);
1439 28 : if (linit_get_type(L) == t_LDESC_PRODUCT)
1440 : {
1441 21 : GEN v = lfunprod_get_fact(linit_get_tech(L)), F = gel(v,1);
1442 21 : long i, l = lg(F);
1443 21 : C = cgetg(l, t_VEC);
1444 70 : for (i = 1; i < l; ++i)
1445 49 : gel(C, i) = zv_to_ZV( lfuncost(gel(F,i), dom, der, bitprec) );
1446 21 : return gc_upto(av, C);
1447 : }
1448 : }
1449 28 : if (!dom) pari_err_TYPE("lfuncost [missing s domain]", L);
1450 28 : C = lfuncost(L,dom,der,bitprec);
1451 28 : return gc_upto(av, zv_to_ZV(C));
1452 : }
1453 :
1454 : static int
1455 10185 : is_dirichlet(GEN ldata)
1456 : {
1457 10185 : switch(ldata_get_type(ldata))
1458 : {
1459 1330 : case t_LFUN_ZETA:
1460 : case t_LFUN_KRONECKER:
1461 1330 : case t_LFUN_CHIZ: return 1;
1462 980 : case t_LFUN_CHIGEN: return ldata_get_degree(ldata)==1;
1463 7875 : default: return 0;
1464 : }
1465 : }
1466 :
1467 : static ulong
1468 11016 : lfuninit_cutoff(GEN ldata)
1469 : {
1470 11016 : GEN gN = ldata_get_conductor(ldata);
1471 : ulong L, N;
1472 11016 : if (ldata_get_type(ldata) == t_LFUN_NF) /* N ~ f^(d-1), exact for d prime */
1473 742 : gN = sqrtnint(gN, ldata_get_degree(ldata) - 1);
1474 11016 : N = itou_or_0(gN);
1475 11016 : if (N > 1000) L = 7000;
1476 11002 : else if (N > 100) L = 5000;
1477 8013 : else if (N > 15) L = 3000;
1478 7558 : else L = 2500;
1479 11016 : return L;
1480 : }
1481 :
1482 : GEN
1483 37349 : lfuninit(GEN lmisc, GEN dom, long der, long bitprec)
1484 : {
1485 37349 : pari_sp av = avma;
1486 : GEN poqk, AB, R, h, theta, ldata, eno, r, domain, tech, k;
1487 : struct lfunp S;
1488 :
1489 37349 : if (is_linit(lmisc))
1490 : {
1491 25316 : long t = linit_get_type(lmisc);
1492 25316 : if (t == t_LDESC_INIT || t == t_LDESC_PRODUCT)
1493 : {
1494 25267 : if (checklfuninit(lmisc, dom, der, bitprec)) return lmisc;
1495 7 : pari_warn(warner,"lfuninit: insufficient initialization");
1496 : }
1497 : }
1498 12089 : ldata = lfunmisc_to_ldata_shallow(lmisc);
1499 :
1500 12089 : switch (ldata_get_type(ldata))
1501 : {
1502 630 : case t_LFUN_NF:
1503 : {
1504 630 : GEN T = gel(ldata_get_an(ldata), 2);
1505 630 : return gc_GEN(av, lfunzetakinit(T, dom, der, bitprec));
1506 : }
1507 91 : case t_LFUN_ABELREL:
1508 : {
1509 91 : GEN T = gel(ldata_get_an(ldata), 2);
1510 91 : return gc_GEN(av, lfunabelianrelinit(gel(T,1), gel(T,2), dom, der, bitprec));
1511 : }
1512 : }
1513 11368 : k = ldata_get_k(ldata);
1514 11368 : parse_dom(gtodouble(k), dom, &S);
1515 : /* Reduce domain for Dirichlet characters. NOT for Abelian t_LFUN_NF,
1516 : * handled above. */
1517 11368 : if (S.dw >= 0 && (!der && is_dirichlet(ldata)))
1518 1750 : S.dh = mindd(S.dh, lfuninit_cutoff(ldata));
1519 11368 : if (S.dw < 0)
1520 : {
1521 98 : if (der)
1522 0 : pari_err_IMPL("domain = [] for derivatives in lfuninit");
1523 98 : if (!is_dirichlet(ldata))
1524 0 : pari_err_IMPL("domain = [] for L functions of degree > 1");
1525 98 : return gc_GEN(av, lfunnoinit(ldata, bitprec));
1526 : }
1527 :
1528 11270 : lfunp_set(ldata, der, bitprec, &S);
1529 11270 : ldata = ldata_newprec(ldata, nbits2prec(S.Dmax));
1530 11270 : r = ldata_get_residue(ldata);
1531 11270 : k = ldata_get_k(ldata); /* if k is t_REAL, ldata_newprec may change it */
1532 : /* Note: all guesses should already have been performed (thetainit more
1533 : * expensive than needed: should be either tdom = 1 or bitprec = S.D).
1534 : * BUT if the root number / polar part do not have an algebraic
1535 : * expression, there is no way to do this until we know the
1536 : * precision, i.e. now. So we can't remove guessing code from here and
1537 : * lfun_init_theta */
1538 11270 : if (r && isintzero(r)) eno = NULL;
1539 : else
1540 : {
1541 11270 : eno = ldata_get_rootno(ldata);
1542 11270 : if (isintzero(eno)) eno = NULL;
1543 : }
1544 11270 : theta = lfun_init_theta(ldata, eno, &S);
1545 11256 : if (eno && !r)
1546 4599 : R = gen_0;
1547 : else
1548 : {
1549 6657 : GEN v = lfunrootres(theta, S.D);
1550 6657 : ldata = shallowcopy(ldata);
1551 6657 : gel(ldata, 6) = gel(v,3);
1552 6657 : r = gel(v,1);
1553 6657 : R = gel(v,2);
1554 6657 : if (isintzero(r)) setlg(ldata,7); else gel(ldata, 7) = r;
1555 : }
1556 11256 : h = divru(mplog2(S.precmax), S.m0);
1557 : /* exp(kh/2 . [0..M]) */
1558 11256 : poqk = gequal0(k) ? NULL
1559 11256 : : gpowers(gprec_w(mpexp(gmul2n(gmul(k,h), -1)), S.precmax), S.M);
1560 11256 : AB = lfuninit_ab(theta, h, &S);
1561 11256 : if (S.bn)
1562 : {
1563 154 : GEN A = gel(AB,1), B = gel(AB,2);
1564 154 : A = lfuninit_pol(A, poqk, S.precmax);
1565 154 : B = lfuninit_pol(B, poqk, S.precmax);
1566 154 : AB = mkvec2(A, B);
1567 : }
1568 : else
1569 11102 : AB = lfuninit_pol(AB, poqk, S.precmax);
1570 11256 : tech = mkvec3(h, AB, R);
1571 11256 : domain = mkvec2(dom, mkvecsmall2(der, bitprec));
1572 11256 : return gc_GEN(av, lfuninit_make(t_LDESC_INIT, ldata, tech, domain));
1573 : }
1574 :
1575 : GEN
1576 553 : lfuninit0(GEN lmisc, GEN dom, long der, long bitprec)
1577 : {
1578 553 : GEN z = lfuninit(lmisc, dom, der, bitprec);
1579 553 : return z == lmisc? gcopy(z): z;
1580 : }
1581 :
1582 : /* If s is a pole of Lambda, return polar part at s; else return NULL */
1583 : static GEN
1584 5431 : lfunpoleresidue(GEN R, GEN s)
1585 : {
1586 : long j;
1587 15768 : for (j = 1; j < lg(R); j++)
1588 : {
1589 10897 : GEN Rj = gel(R, j), be = gel(Rj, 1);
1590 10897 : if (gequal(s, be)) return gel(Rj, 2);
1591 : }
1592 4871 : return NULL;
1593 : }
1594 :
1595 : /* Compute contribution of polar part at s when not a pole. */
1596 : static GEN
1597 9021 : veccothderivn(GEN a, long n)
1598 : {
1599 : long i;
1600 9021 : pari_sp av = avma;
1601 9021 : GEN c = pol_x(0), cp = mkpoln(3, gen_m1, gen_0, gen_1);
1602 9021 : GEN v = cgetg(n+2, t_VEC);
1603 9021 : gel(v, 1) = poleval(c, a);
1604 27182 : for(i = 2; i <= n+1; i++)
1605 : {
1606 18161 : c = ZX_mul(ZX_deriv(c), cp);
1607 18161 : gel(v, i) = gdiv(poleval(c, a), mpfact(i-1));
1608 : }
1609 9021 : return gc_GEN(av, v);
1610 : }
1611 :
1612 : static GEN
1613 9140 : polepart(long n, GEN h, GEN C)
1614 : {
1615 9140 : GEN h2n = gpowgs(gdiv(h, gen_2), n-1);
1616 9140 : GEN res = gmul(h2n, gel(C,n));
1617 9140 : return odd(n)? res : gneg(res);
1618 : }
1619 :
1620 : static GEN
1621 4381 : lfunsumcoth(GEN R, GEN s, GEN h, long prec)
1622 : {
1623 : long i,j;
1624 4381 : GEN S = gen_0;
1625 13402 : for (j = 1; j < lg(R); ++j)
1626 : {
1627 9021 : GEN r = gel(R,j), be = gel(r,1), Rj = gel(r, 2);
1628 9021 : long e = valser(Rj);
1629 9021 : GEN z1 = gexpm1(gmul(h, gsub(s,be)), prec); /* exp(h(s-beta))-1 */
1630 9021 : GEN c1 = gaddgs(gdivsg(2, z1), 1); /* coth((h/2)(s-beta)) */
1631 9021 : GEN C1 = veccothderivn(c1, 1-e);
1632 18161 : for (i = e; i < 0; i++)
1633 : {
1634 9140 : GEN Rbe = mysercoeff(Rj, i);
1635 9140 : GEN p1 = polepart(-i, h, C1);
1636 9140 : S = gadd(S, gmul(Rbe, p1));
1637 : }
1638 : }
1639 4381 : return gmul2n(S, -1);
1640 : }
1641 :
1642 : static GEN lfunlambda_OK(GEN linit, GEN s, GEN sdom, long bitprec);
1643 : /* L is a t_LDESC_PRODUCT or t_LDESC_INIT Linit */
1644 : static GEN
1645 2456 : _product(GEN (*fun)(GEN,GEN,long), GEN L, GEN s, long bitprec)
1646 : {
1647 2456 : GEN ldata = linit_get_ldata(L), v, r, F, E, C, cs;
1648 : long i, l;
1649 : int isreal;
1650 2456 : if (linit_get_type(L) == t_LDESC_INIT) return fun(ldata, s, bitprec);
1651 1924 : v = lfunprod_get_fact(linit_get_tech(L));
1652 1924 : F = gel(v,1); E = gel(v,2); C = gel(v,3); l = lg(F);
1653 1924 : cs = conj_i(s); isreal = gequal(imag_i(s), imag_i(cs));
1654 6374 : for (i = 1, r = gen_1; i < l; ++i)
1655 : {
1656 4450 : GEN f = fun(gel(F, i), s, bitprec);
1657 4450 : if (typ(f)==t_VEC) f = RgV_prod(f);
1658 4450 : if (E[i]) r = gmul(r, gpowgs(f, E[i]));
1659 4450 : if (C[i])
1660 : {
1661 0 : GEN fc = isreal? f: conj_i(fun(gel(F, i), cs, bitprec));
1662 0 : r = gmul(r, gpowgs(fc, C[i]));
1663 : }
1664 : }
1665 1924 : return (ldata_isreal(ldata) && gequal0(imag_i(s)))? real_i(r): r;
1666 : }
1667 :
1668 : /* s a t_SER; # terms - 1 */
1669 : static long
1670 2318 : der_level(GEN s)
1671 2318 : { return signe(s)? lg(s)-3: valser(s)-1; }
1672 :
1673 : /* s a t_SER; return coeff(s, X^0) */
1674 : static GEN
1675 1260 : ser_coeff0(GEN s) { return simplify_shallow(polcoef_i(s, 0, -1)); }
1676 :
1677 : static GEN
1678 20480 : get_domain(GEN s, GEN *dom, long *der)
1679 : {
1680 20480 : GEN sa = s;
1681 20480 : *der = 0;
1682 20480 : switch(typ(s))
1683 : {
1684 7 : case t_POL:
1685 7 : case t_RFRAC: s = toser_i(s);
1686 1260 : case t_SER:
1687 1260 : *der = der_level(s);
1688 1260 : sa = ser_coeff0(s);
1689 : }
1690 20480 : *dom = mkvec3(real_i(sa), gen_0, gabs(imag_i(sa),DEFAULTPREC));
1691 20480 : return s;
1692 : }
1693 : /* assume s went through get_domain and s/bitprec belong to domain */
1694 : static GEN
1695 34192 : lfunlambda_OK(GEN linit, GEN s, GEN sdom, long bitprec)
1696 : {
1697 34192 : GEN dom, eno, ldata, tech, h, pol, k2, cost, S, S0 = NULL;
1698 : long prec, prec0;
1699 : struct lfunp D, D0;
1700 :
1701 34192 : if (linit_get_type(linit) == t_LDESC_PRODUCT)
1702 1679 : return _product(&lfunlambda, linit, s, bitprec);
1703 32513 : ldata = linit_get_ldata(linit);
1704 32513 : eno = ldata_get_rootno(ldata);
1705 32513 : tech = linit_get_tech(linit);
1706 32513 : dom = lfun_get_dom(tech);
1707 32513 : if (lg(dom) == 1) return lfunlambda(linit, s, bitprec); /* FIXME:not OK! */
1708 32513 : h = lfun_get_step(tech); prec = realprec(h);
1709 : /* try to reduce accuracy */
1710 32513 : parse_dom(0, sdom, &D0);
1711 32513 : parse_dom(0, dom, &D);
1712 32513 : if (0.8 * D.dh > D0.dh)
1713 : {
1714 16270 : cost = lfuncost(linit, sdom, typ(s)==t_SER? der_level(s): 0, bitprec);
1715 16270 : prec0 = nbits2prec(cost[2]);
1716 16270 : if (prec0 < prec) { prec = prec0; h = gprec_w(h, prec); }
1717 : }
1718 32513 : pol = lfun_get_pol(tech);
1719 32513 : s = gprec_w(s, prec);
1720 32513 : if (ldata_get_residue(ldata))
1721 : {
1722 4780 : GEN R = lfun_get_Residue(tech);
1723 4780 : GEN Ra = lfunpoleresidue(R, s);
1724 4780 : if (Ra) return gprec_w(Ra, nbits2prec(bitprec));
1725 4381 : S0 = lfunsumcoth(R, s, h, prec);
1726 : }
1727 32114 : k2 = lfun_get_k2(tech);
1728 32114 : if (typ(pol)==t_POL && typ(s) != t_SER && gequal(real_i(s), k2))
1729 25485 : { /* on critical line: shortcut */
1730 25485 : GEN polz, b = imag_i(s);
1731 25485 : polz = gequal0(b)? poleval(pol,gen_1): poleval(pol, expIr(gmul(h,b)));
1732 25485 : S = gadd(polz, gmulvec(eno, conj_i(polz)));
1733 : }
1734 : else
1735 : {
1736 6629 : GEN z = gexp(gmul(h, gsub(s, k2)), prec);
1737 6629 : GEN zi = ginv(z), zc = conj_i(zi);
1738 6629 : if (typ(pol)==t_POL)
1739 6433 : S = gadd(poleval(pol, z), gmulvec(eno, conj_i(poleval(pol, zc))));
1740 : else
1741 196 : S = gadd(poleval(gel(pol,1), z), gmulvec(eno, poleval(gel(pol,2), zi)));
1742 : }
1743 32114 : if (S0) S = gadd(S,S0);
1744 32114 : return gprec_w(gmul(S,h), nbits2prec(bitprec));
1745 : }
1746 :
1747 : static long
1748 38973 : lfunspec_OK(GEN lmisc, GEN s, GEN *pldata)
1749 : {
1750 38973 : long t, large = 0;
1751 : GEN ldata;
1752 38973 : *pldata = ldata = lfunmisc_to_ldata_shallow(lmisc);
1753 38966 : if (!is_linit(lmisc)) lmisc = ldata;
1754 25981 : else switch(linit_get_type(lmisc))
1755 : {
1756 25939 : case t_LDESC_INIT: case t_LDESC_PRODUCT:
1757 25939 : if (lg(lfun_get_dom(linit_get_tech(lmisc))) == 1) large = 1;
1758 25939 : break;
1759 42 : default: return 0;
1760 : }
1761 38924 : switch(typ(s))
1762 : {
1763 38553 : case t_INT: case t_REAL: case t_FRAC: case t_COMPLEX: break;
1764 371 : default: return 0;
1765 : }
1766 38553 : t = ldata_get_type(ldata);
1767 38553 : switch(t)
1768 : {
1769 7930 : case t_LFUN_KRONECKER: case t_LFUN_ZETA:
1770 7930 : if (typ(s) == t_INT && !is_bigint(s)) return 1;
1771 : /* fall through */
1772 : case t_LFUN_NF: case t_LFUN_CHIZ:
1773 8056 : if (!large)
1774 7762 : large = (fabs(gtodouble(imag_i(s))) >= lfuninit_cutoff(ldata));
1775 8056 : break;
1776 2057 : case t_LFUN_CHIGEN:
1777 2057 : if (ldata_get_degree(ldata) != 1) return 0;
1778 1735 : if (!large)
1779 1483 : large = (fabs(gtodouble(imag_i(s))) >= lfuninit_cutoff(ldata));
1780 1735 : break;
1781 : }
1782 34150 : if (large)
1783 : {
1784 1904 : if (t == t_LFUN_NF)
1785 : {
1786 14 : GEN an = ldata_get_an(ldata), nf = gel(an,2), G = galoisinit(nf, NULL);
1787 14 : if (isintzero(G) || !group_isabelian(galois_group(G))) return 0;
1788 : }
1789 1904 : return 2;
1790 : }
1791 32246 : return 0;
1792 : }
1793 :
1794 : GEN
1795 5136 : lfunlambda(GEN lmisc, GEN s, long bitprec)
1796 : {
1797 5136 : pari_sp av = avma;
1798 5136 : GEN linit = NULL, dom, z;
1799 : long der;
1800 5136 : s = get_domain(s, &dom, &der);
1801 5136 : if (!der)
1802 : {
1803 : GEN ldata;
1804 4289 : long t = lfunspec_OK(lmisc, s, &ldata);
1805 4289 : if (t == 1)
1806 : { /* special value ? */
1807 553 : GEN z = lfun(ldata, s, bitprec), gv = ldata_get_gammavec(ldata);
1808 553 : long e = itou(gel(gv, 1));
1809 553 : if (!isintzero(z) && (e || gsigne(s) > 0)) /* TODO */
1810 : {
1811 476 : GEN q = ldata_get_conductor(ldata);
1812 476 : long prec = nbits2prec(bitprec);
1813 476 : GEN se, r, Q = divir(q, mppi(prec));
1814 476 : se = gmul2n(gaddgs(s, e), -1);
1815 476 : r = gmul(gpow(Q, se, prec), ggamma(se, prec));
1816 476 : if (e && !equali1(q)) r = gdiv(r, gsqrt(q, prec));
1817 504 : return gc_upto(av, gmul(r, z));
1818 : }
1819 : }
1820 3813 : if (is_linit(lmisc)) linit = lmisc; else lmisc = ldata;
1821 3813 : if (t == 2)
1822 28 : return gc_GEN(av, linit? _product(&lfunlambda, linit, s, bitprec)
1823 0 : : lfunlambdalarge(ldata, s, bitprec));
1824 : }
1825 4632 : linit = lfuninit(lmisc, dom, der, bitprec);
1826 4632 : z = lfunlambda_OK(linit,s, dom, bitprec);
1827 4632 : return gc_GEN(av, z);
1828 : }
1829 :
1830 : static long
1831 19481 : is_ser(GEN x)
1832 : {
1833 19481 : long t = typ(x);
1834 19481 : if (t == t_SER) return 1;
1835 17430 : if (!is_vec_t(t) || lg(x)==1) return 0;
1836 350 : if (typ(gel(x,1))==t_SER) return 1;
1837 252 : return 0;
1838 : }
1839 :
1840 : static GEN
1841 371 : lfunser(GEN L)
1842 : {
1843 371 : long v = valser(L);
1844 371 : if (v > 0) return gen_0;
1845 329 : if (v == 0) L = gel(L, 2);
1846 : else
1847 203 : setlg(L, minss(lg(L), 2-v));
1848 329 : return L;
1849 : }
1850 :
1851 : static GEN
1852 371 : lfunservec(GEN x)
1853 : {
1854 371 : if (typ(x)==t_SER) return lfunser(x);
1855 0 : pari_APPLY_same(lfunser(gel(x,i)))
1856 : }
1857 : static GEN
1858 105 : lfununext(GEN L)
1859 : {
1860 105 : setlg(L, maxss(lg(L)-1, valser(L)? 2: 3));
1861 105 : return normalizeser(L);
1862 : }
1863 : static GEN
1864 105 : lfununextvec(GEN x)
1865 : {
1866 105 : if (typ(x)==t_SER) return lfununext(x);
1867 0 : pari_APPLY_same(lfununext(gel(x,i)));
1868 : }
1869 :
1870 : /* assume lmisc is an linit, s went through get_domain and s/bitprec belong
1871 : * to domain */
1872 : static GEN
1873 10402 : lfun_OK(GEN linit, GEN s, GEN sdom, long bitprec)
1874 : {
1875 10402 : GEN N, gas, S, FVga, res, ss = s;
1876 10402 : long prec = nbits2prec(bitprec), ext;
1877 :
1878 10402 : FVga = lfun_get_factgammavec(linit_get_tech(linit));
1879 10402 : S = lfunlambda_OK(linit, s, sdom, bitprec);
1880 10402 : if (is_ser(S))
1881 : {
1882 1715 : GEN r = typ(S)==t_SER ? S : gel(S,1);
1883 1715 : long d = lg(r) - 2 + fracgammadegree(FVga);
1884 1715 : if (typ(s) == t_SER)
1885 1386 : ss = sertoser(s, d);
1886 : else
1887 329 : ss = deg1ser_shallow(gen_1, s, varn(r), d);
1888 : }
1889 10402 : gas = gammafactproduct(FVga, ss, &ext, prec);
1890 10402 : N = ldata_get_conductor(linit_get_ldata(linit));
1891 10402 : res = gdiv(S, gmul(gpow(N, gdivgu(ss, 2), prec), gas));
1892 10402 : if (typ(s) != t_SER && is_ser(res)) res = lfunservec(res);
1893 10031 : else if (ext) res = lfununextvec(res);
1894 10402 : return gprec_w(res, prec);
1895 : }
1896 :
1897 : GEN
1898 13643 : lfun(GEN lmisc, GEN s, long bitprec)
1899 : {
1900 13643 : pari_sp av = avma;
1901 13643 : GEN linit = NULL, ldata, dom, z;
1902 : long der;
1903 13643 : s = get_domain(s, &dom, &der);
1904 13643 : if (der && typ(s) != t_SER)
1905 : {
1906 0 : if (lfunspec_OK(lmisc, s, &ldata))
1907 : {
1908 0 : linit = lfuninit(lmisc, cgetg(1,t_VEC), 0, bitprec);
1909 0 : return derivnumk((void*)linit, (GEN(*)(void*,GEN,long))&lfun,
1910 : s, stoi(der), nbits2prec(bitprec));
1911 : }
1912 : }
1913 : else
1914 : {
1915 13643 : long t = lfunspec_OK(lmisc, s, &ldata);
1916 13636 : if (t == 1)
1917 : { /* special value ? */
1918 3353 : long D = itos_or_0(gel(ldata_get_an(ldata), 2)), ss = itos(s);
1919 3353 : if (D)
1920 : {
1921 3353 : if (ss <= 0) return lfunquadneg(D, ss);
1922 : /* ss > 0 */
1923 770 : if ((!odd(ss) && D > 0) || (odd(ss) && D < 0))
1924 : {
1925 672 : long prec = nbits2prec(bitprec), q = labs(D);
1926 672 : ss = 1 - ss; /* <= 0 */
1927 672 : z = powrs(divrs(mppi(prec + EXTRAPREC64), q), 1-ss);
1928 672 : z = mulrr(shiftr(z, -ss), sqrtr_abs(utor(q, prec)));
1929 672 : z = gdiv(z, mpfactr(-ss, prec));
1930 672 : if (smodss(ss, 4) > 1) togglesign(z);
1931 672 : return gmul(z, lfunquadneg(D, ss));
1932 : }
1933 : }
1934 : }
1935 10381 : if (is_linit(lmisc)) linit = lmisc; else lmisc = ldata;
1936 10381 : if (t == 2)
1937 1211 : return gc_GEN(av, linit? _product(&lfun, linit, s, bitprec)
1938 231 : : lfunlarge(ldata, s, bitprec));
1939 : }
1940 9401 : linit = lfuninit(lmisc, dom, der, bitprec);
1941 9387 : z = lfun_OK(linit, s, dom, bitprec);
1942 9387 : return gc_GEN(av, z);
1943 : }
1944 :
1945 : /* given a t_SER a+x*s(x), return x*s(x), shallow */
1946 : static GEN
1947 42 : sersplit1(GEN s, GEN *head)
1948 : {
1949 42 : long i, l = lg(s);
1950 : GEN y;
1951 42 : *head = simplify_shallow(mysercoeff(s, 0));
1952 42 : if (valser(s) > 0) return s;
1953 28 : y = cgetg(l-1, t_SER); y[1] = s[1];
1954 28 : setvalser(y, 1);
1955 140 : for (i=3; i < l; i++) gel(y,i-1) = gel(s,i);
1956 28 : return normalizeser(y);
1957 : }
1958 :
1959 : /* order of pole of Lambda at s (0 if regular point) */
1960 : static long
1961 2310 : lfunlambdaord(GEN linit, GEN s)
1962 : {
1963 2310 : GEN tech = linit_get_tech(linit);
1964 2310 : if (linit_get_type(linit)==t_LDESC_PRODUCT)
1965 : {
1966 287 : GEN v = lfunprod_get_fact(linit_get_tech(linit));
1967 287 : GEN F = gel(v, 1), E = gel(v, 2), C = gel(v, 3);
1968 287 : long i, ex = 0, l = lg(F);
1969 980 : for (i = 1; i < l; i++)
1970 693 : ex += lfunlambdaord(gel(F,i), s) * (E[i]+C[i]);
1971 287 : return ex;
1972 : }
1973 2023 : if (ldata_get_residue(linit_get_ldata(linit)))
1974 : {
1975 651 : GEN r = lfunpoleresidue(lfun_get_Residue(tech), s);
1976 651 : if (r) return lg(r)-2;
1977 : }
1978 1862 : return 0;
1979 : }
1980 :
1981 : static GEN
1982 126 : derser(GEN res, long m)
1983 : {
1984 126 : long v = valser(res);
1985 126 : if (v > m) return gen_0;
1986 126 : if (v >= 0)
1987 126 : return gmul(mysercoeff(res, m), mpfact(m));
1988 : else
1989 0 : return derivn(res, m, -1);
1990 : }
1991 :
1992 : static GEN
1993 189 : derservec(GEN x, long m) { pari_APPLY_same(derser(gel(x,i),m)) }
1994 :
1995 : /* derivative of order m > 0 of L (flag = 0) or Lambda (flag = 1) */
1996 : static GEN
1997 1708 : lfunderiv(GEN lmisc, long m, GEN s, long flag, long bitprec)
1998 : {
1999 1708 : pari_sp ltop = avma;
2000 1708 : GEN res, S = NULL, linit, ldata, dom;
2001 1708 : long der, prec = nbits2prec(bitprec);
2002 1708 : if (m <= 0) pari_err_DOMAIN("lfun", "D", "<=", gen_0, stoi(m));
2003 1701 : s = get_domain(s, &dom, &der);
2004 1701 : if (typ(s) != t_SER && lfunspec_OK(lmisc, s, &ldata) == 2)
2005 : {
2006 28 : linit = lfuninit(lmisc, cgetg(1,t_VEC), 0, bitprec);
2007 28 : return derivnumk((void*)linit, (GEN(*)(void*,GEN,long))&lfun,
2008 : s, stoi(der + m), prec);
2009 : }
2010 1673 : linit = lfuninit(lmisc, dom, der + m, bitprec);
2011 1673 : if (lg(lfun_get_dom(linit_get_tech(linit))) == 1)
2012 14 : pari_err_IMPL("domain = [] for derivatives in lfuninit");
2013 1659 : if (typ(s) == t_SER)
2014 : {
2015 : GEN a;
2016 42 : if (valser(s) < 0) pari_err_DOMAIN("lfun","valuation", "<", gen_0, s);
2017 42 : S = sersplit1(s, &a);
2018 42 : s = deg1ser_shallow(gen_1, a, varn(S), m + ceildivuu(lg(s)-2, valser(S)));
2019 : }
2020 : else
2021 : {
2022 1617 : long e = lfunlambdaord(linit, s) + m + 1;
2023 : /* HACK: pretend lfuninit was done to right accuracy */
2024 1617 : if (gequal0(s)) { s = gen_0; e--; }
2025 1617 : s = deg1ser_shallow(gen_1, s, 0, e);
2026 : }
2027 1659 : res = flag ? lfunlambda_OK(linit, s, dom, bitprec):
2028 1015 : lfun_OK(linit, s, dom, bitprec);
2029 1659 : if (S)
2030 42 : res = gsubst(derivn(res, m, -1), varn(S), S);
2031 1617 : else if (typ(res)==t_SER)
2032 : {
2033 1554 : long v = valser(res);
2034 1554 : if (v > m) { set_avma(ltop); return gen_0; }
2035 1540 : if (v >= 0)
2036 1414 : res = gmul(mysercoeff(res, m), mpfact(m));
2037 : else
2038 126 : res = derivn(res, m, -1);
2039 : }
2040 63 : else if (is_ser(res))
2041 63 : res = derservec(res, m);
2042 1645 : return gc_GEN(ltop, gprec_w(res, prec));
2043 : }
2044 :
2045 : GEN
2046 1596 : lfunlambda0(GEN lmisc, GEN s, long der, long bitprec)
2047 : {
2048 665 : return der? lfunderiv(lmisc, der, s, 1, bitprec)
2049 2254 : : lfunlambda(lmisc, s, bitprec);
2050 : }
2051 :
2052 : GEN
2053 7392 : lfun0(GEN lmisc, GEN s, long der, long bitprec)
2054 : {
2055 1043 : return der? lfunderiv(lmisc, der, s, 0, bitprec)
2056 8421 : : lfun(lmisc, s, bitprec);
2057 : }
2058 :
2059 : GEN
2060 19389 : lfunhardy(GEN lmisc, GEN t, long bitprec)
2061 : {
2062 19389 : pari_sp ltop = avma;
2063 19389 : long prec = nbits2prec(bitprec), d, isbig = 0;
2064 : GEN linit, h, ldata, tech, w2, k2, E, a, argz, z;
2065 :
2066 19389 : switch(typ(t))
2067 : {
2068 19382 : case t_INT: case t_FRAC: case t_REAL: break;
2069 7 : default: pari_err_TYPE("lfunhardy",t);
2070 : }
2071 19382 : if (lfunspec_OK(lmisc, mkcomplex(gen_0, t), &ldata) == 2)
2072 : {
2073 868 : long B = bitprec + maxss(gexpo(t), 0);
2074 868 : GEN L = NULL;
2075 868 : isbig = 1;
2076 868 : k2 = ghalf;
2077 868 : z = mkcomplex(k2, t);
2078 868 : if (is_linit(lmisc))
2079 : {
2080 742 : linit = lmisc;
2081 742 : if (linit_get_type(linit) == t_LDESC_PRODUCT)
2082 14 : L = mkvec(linit);/*HACK*/
2083 : }
2084 : else
2085 : {
2086 126 : linit = lfunnoinit(ldata, B);
2087 126 : ldata = linit_get_ldata(linit); /* make sure eno is included */
2088 : }
2089 868 : h = lfunloglambdalarge(L? L: ldata, gprec_w(z, nbits2prec(B)), B);
2090 868 : tech = linit_get_tech(linit);
2091 : }
2092 : else
2093 : {
2094 18514 : GEN k = ldata_get_k(ldata);
2095 18514 : GEN dom = mkvec3(gmul2n(k, -1), gen_0, gabs(t,LOWDEFAULTPREC));
2096 18514 : if (!is_linit(lmisc)) lmisc = ldata;
2097 18514 : linit = lfuninit(lmisc, dom, 0, bitprec);
2098 18514 : tech = linit_get_tech(linit);
2099 18514 : k2 = lfun_get_k2(tech);
2100 18514 : z = mkcomplex(k2, t);
2101 18514 : h = lfunlambda_OK(linit, z, dom, bitprec);
2102 : }
2103 19382 : w2 = lfun_get_w2(tech);
2104 19382 : E = lfun_get_expot(tech); /* 4E = d(k2 - 1) + real(vecsum(Vga)) */
2105 19382 : d = ldata_get_degree(ldata);
2106 : /* more accurate than garg: k/2 in Q */
2107 19382 : argz = gequal0(k2)? Pi2n(-1, prec): gatan(gdiv(t, k2), prec);
2108 19382 : prec = precision(argz);
2109 : /* prec may have increased: don't lose accuracy if |z|^2 is exact */
2110 19382 : a = gsub(gmulsg(d, gmul(t, gmul2n(argz,-1))),
2111 : gmul(E, glog(gnorm(z),prec)));
2112 19382 : if (!isint1(w2) && typ(ldata_get_dual(ldata))==t_INT)
2113 16198 : h = isbig ? gadd(h, glog(w2, prec)) : mulrealvec(h, w2);
2114 19382 : if (typ(h) == t_COMPLEX && gexpo(imag_i(h)) < -(bitprec >> 1))
2115 2575 : h = real_i(h);
2116 19382 : if (isbig) h = greal(gexp(gadd(h, a), prec));
2117 18514 : else h = gmul(h, gexp(a, prec));
2118 19382 : return gc_upto(ltop, h);
2119 : }
2120 :
2121 : /* L = log(t); return \sum_{i = 0}^{v-1} R[-i-1] L^i/i! */
2122 : static GEN
2123 2086 : theta_pole_contrib(GEN R, long v, GEN L)
2124 : {
2125 2086 : GEN s = mysercoeff(R,-v);
2126 : long i;
2127 2191 : for (i = v-1; i >= 1; i--)
2128 105 : s = gadd(mysercoeff(R,-i), gdivgu(gmul(s,L), i));
2129 2086 : return s;
2130 : }
2131 : /* subtract successively rather than adding everything then subtracting.
2132 : * The polar part is "large" and suffers from cancellation: a little stabler
2133 : * this way */
2134 : static GEN
2135 7455 : theta_add_polar_part(GEN S, GEN R, GEN t, long prec)
2136 : {
2137 7455 : GEN logt = NULL;
2138 7455 : long j, l = lg(R);
2139 9541 : for (j = 1; j < l; j++)
2140 : {
2141 2086 : GEN Rj = gel(R,j), b = gel(Rj,1), Rb = gel(Rj,2);
2142 2086 : long v = -valser(Rb);
2143 2086 : if (v > 1 && !logt) logt = glog(t, prec);
2144 2086 : S = gsub(S, gmul(theta_pole_contrib(Rb,v,logt), gpow(t,b,prec)));
2145 : }
2146 7455 : return S;
2147 : }
2148 :
2149 : static long
2150 3689 : lfuncheckfeq_i(GEN theta, GEN thetad, GEN t0, GEN t0i, long bitprec)
2151 : {
2152 3689 : GEN ldata = linit_get_ldata(theta);
2153 : GEN S0, S0i, w, eno;
2154 3689 : long prec = nbits2prec(bitprec);
2155 3689 : if (thetad)
2156 70 : S0 = lfuntheta(thetad, t0, 0, bitprec);
2157 : else
2158 3619 : S0 = conj_i(lfuntheta(theta, conj_i(t0), 0, bitprec));
2159 3689 : S0i = lfuntheta(theta, t0i, 0, bitprec);
2160 :
2161 3689 : eno = ldata_get_rootno(ldata);
2162 3689 : if (ldata_get_residue(ldata))
2163 : {
2164 1015 : GEN R = theta_get_R(linit_get_tech(theta));
2165 1015 : if (gequal0(R))
2166 : {
2167 : GEN v, r;
2168 105 : long t = ldata_get_type(ldata);
2169 105 : if (t == t_LFUN_NF || t == t_LFUN_ABELREL)
2170 : { /* inefficient since theta not needed; no need to optimize for this
2171 : (artificial) query [e.g. lfuncheckfeq(t_POL)] */
2172 42 : GEN L = lfuninit(ldata,zerovec(3),0,bitprec);
2173 42 : return lfuncheckfeq(L,t0,bitprec);
2174 : }
2175 63 : v = lfunrootres(theta, bitprec);
2176 63 : r = gel(v,1);
2177 63 : if (gequal0(eno)) eno = gel(v,3);
2178 63 : R = lfunrtoR_i(ldata, r, eno, nbits2prec(bitprec));
2179 : }
2180 973 : S0i = theta_add_polar_part(S0i, R, t0, prec);
2181 : }
2182 3647 : if (gequal0(S0i) || gequal0(S0)) pari_err_PREC("lfuncheckfeq");
2183 :
2184 3647 : w = gdivvec(S0i, gmul(S0, gpow(t0, ldata_get_k(ldata), prec)));
2185 : /* missing rootno: guess it */
2186 3647 : if (gequal0(eno)) eno = lfunrootno(theta, bitprec);
2187 3647 : w = gsubvec(w, eno);
2188 3647 : if (thetad) w = gdivvec(w, eno); /* |eno| may be large in non-dual case */
2189 3647 : return gexpo(w);
2190 : }
2191 :
2192 : /* Check whether the coefficients, conductor, weight, polar part and root
2193 : * number are compatible with the functional equation at t0 and 1/t0.
2194 : * Different from lfunrootres. */
2195 : long
2196 3822 : lfuncheckfeq(GEN lmisc, GEN t0, long bitprec)
2197 : {
2198 : GEN ldata, theta, thetad, t0i;
2199 : pari_sp av;
2200 :
2201 3822 : if (is_linit(lmisc) && linit_get_type(lmisc)==t_LDESC_PRODUCT)
2202 : {
2203 168 : GEN v = lfunprod_get_fact(linit_get_tech(lmisc)), F = gel(v,1);
2204 168 : long i, b = -bitprec, l = lg(F);
2205 560 : for (i = 1; i < l; i++) b = maxss(b, lfuncheckfeq(gel(F,i), t0, bitprec));
2206 168 : return b;
2207 : }
2208 3654 : av = avma;
2209 3654 : if (!t0)
2210 : { /* ~Pi/3 + I/7, some random complex number */
2211 3479 : t0 = mkcomplex(uutoQ(355,339), uutoQ(1,7));
2212 3479 : t0i = ginv(t0);
2213 : }
2214 175 : else if (gcmpgs(gnorm(t0), 1) < 0) { t0i = t0; t0 = ginv(t0); }
2215 119 : else t0i = ginv(t0);
2216 : /* |t0| >= 1 */
2217 3654 : theta = lfunthetacheckinit(lmisc, t0i, 0, bitprec);
2218 3647 : ldata = linit_get_ldata(theta);
2219 3647 : thetad = theta_dual(theta, ldata_get_dual(ldata));
2220 3647 : return gc_long(av, lfuncheckfeq_i(theta, thetad, t0, t0i, bitprec));
2221 : }
2222 :
2223 : /*******************************************************************/
2224 : /* Compute root number and residues */
2225 : /*******************************************************************/
2226 : /* round root number to \pm 1 if close to integer. */
2227 : static GEN
2228 6825 : ropm1(GEN w, long prec)
2229 : {
2230 : long e;
2231 : GEN r;
2232 6825 : if (typ(w) == t_INT) return w;
2233 6419 : r = grndtoi(w, &e);
2234 6419 : return (e < -prec/2)? r: w;
2235 : }
2236 :
2237 : /* theta for t=1/sqrt(2) and t2==2t simultaneously, saving 25% of the work.
2238 : * Assume correct initialization (no thetacheck) */
2239 : static void
2240 441 : lfunthetaspec(GEN linit, long bitprec, GEN *pv, GEN *pv2)
2241 : {
2242 441 : pari_sp av = avma, av2;
2243 : GEN t, Vga, an, K, ldata, thetainit, v, v2, vroots;
2244 : long L, prec, n, d;
2245 :
2246 441 : ldata = linit_get_ldata(linit);
2247 441 : thetainit = linit_get_tech(linit);
2248 441 : prec = nbits2prec(bitprec);
2249 441 : Vga = ldata_get_gammavec(ldata); d = lg(Vga)-1;
2250 441 : if (Vgaeasytheta(Vga))
2251 : {
2252 224 : GEN v2 = sqrtr(real2n(1, nbits2prec(bitprec)));
2253 224 : GEN v = shiftr(v2,-1);
2254 224 : *pv = lfuntheta(linit, v, 0, bitprec);
2255 224 : *pv2= lfuntheta(linit, v2, 0, bitprec);
2256 224 : return;
2257 : }
2258 217 : an = RgV_kill0( theta_get_an(thetainit) );
2259 217 : L = lg(an)-1;
2260 : /* to compute theta(1/sqrt(2)) */
2261 217 : t = ginv(gsqrt(gmul2n(ldata_get_conductor(ldata), 1), prec));
2262 : /* t = 1/sqrt(2N) */
2263 :
2264 : /* From then on, the code is generic and could be used to compute
2265 : * theta(t) / theta(2t) without assuming t = 1/sqrt(2) */
2266 217 : K = theta_get_K(thetainit);
2267 217 : vroots = mkvroots(d, L, prec);
2268 217 : t = gpow(t, gdivgu(gen_2, d), prec); /* rt variant: t->t^(2/d) */
2269 : /* v = \sum_{n <= L, n odd} a_n K(nt) */
2270 1822800 : for (v = gen_0, n = 1; n <= L; n+=2)
2271 : {
2272 1822583 : GEN tn, Kn, a = gel(an, n);
2273 :
2274 1822583 : if (!a) continue;
2275 119931 : av2 = avma;
2276 119931 : tn = gmul(t, gel(vroots,n));
2277 119931 : Kn = gammamellininvrt(K, tn, bitprec);
2278 119931 : v = gc_upto(av2, gadd(v, gmul(a,Kn)));
2279 : }
2280 : /* v += \sum_{n <= L, n even} a_n K(nt), v2 = \sum_{n <= L/2} a_n K(2n t) */
2281 1822695 : for (v2 = gen_0, n = 1; n <= L/2; n++)
2282 : {
2283 1822478 : GEN t2n, K2n, a = gel(an, n), a2 = gel(an,2*n);
2284 :
2285 1822478 : if (!a && !a2) continue;
2286 126714 : av2 = avma;
2287 126714 : t2n = gmul(t, gel(vroots,2*n));
2288 126714 : K2n = gc_upto(av2, gammamellininvrt(K, t2n, bitprec));
2289 126714 : if (a) v2 = gadd(v2, gmul(a, K2n));
2290 126714 : if (a2) v = gadd(v, gmul(a2,K2n));
2291 : }
2292 217 : *pv = v;
2293 217 : *pv2 = v2;
2294 217 : (void)gc_all(av, 2, pv,pv2);
2295 : }
2296 :
2297 : static GEN
2298 413 : Rtor(GEN a, GEN R, GEN ldata, long prec)
2299 : {
2300 413 : GEN FVga = gammafactor(ldata_get_gammavec(ldata));
2301 413 : GEN Na = gpow(ldata_get_conductor(ldata), gdivgu(a,2), prec);
2302 : long ext;
2303 413 : return gdiv(R, gmul(Na, gammafactproduct(FVga, a, &ext, prec)));
2304 : }
2305 :
2306 : /* v = theta~(t), vi = theta(1/t) */
2307 : static GEN
2308 6482 : get_eno(GEN R, GEN k, GEN t, GEN v, GEN vi, long vx, long bitprec, long force)
2309 : {
2310 6482 : long prec = nbits2prec(bitprec);
2311 6482 : GEN a0, a1, S = deg1pol_shallow(gmul(gpow(t,k,prec), gneg(v)), vi, vx);
2312 :
2313 6482 : S = theta_add_polar_part(S, R, t, prec);
2314 6482 : if (typ(S) != t_POL || degpol(S) != 1) return NULL;
2315 6482 : a1 = gel(S,3); if (!force && gexpo(a1) < -bitprec/4) return NULL;
2316 6412 : a0 = gel(S,2);
2317 6412 : return gdivvec(a0, gneg(a1));
2318 :
2319 : }
2320 : /* Return w using theta(1/t) - w t^k \bar{theta}(t) = polar_part(t,w).
2321 : * The full Taylor expansion of L must be known */
2322 : GEN
2323 6412 : lfunrootno(GEN linit, long bitprec)
2324 : {
2325 : GEN ldata, t, eno, v, vi, R, thetad;
2326 6412 : long c = 0, prec = nbits2prec(bitprec), vx = fetch_var();
2327 : GEN k;
2328 : pari_sp av;
2329 :
2330 : /* initialize for t > 1/sqrt(2) */
2331 6412 : linit = lfunthetacheckinit(linit, dbltor(sqrt(0.5)), 0, bitprec);
2332 6412 : ldata = linit_get_ldata(linit);
2333 6412 : k = ldata_get_k(ldata);
2334 6426 : R = ldata_get_residue(ldata)? lfunrtoR_eno(ldata, pol_x(vx), prec)
2335 6412 : : cgetg(1, t_VEC);
2336 6412 : t = gen_1;
2337 6412 : v = lfuntheta(linit, t, 0, bitprec);
2338 6412 : thetad = theta_dual(linit, ldata_get_dual(ldata));
2339 6412 : vi = !thetad ? conj_i(v): lfuntheta(thetad, t, 0, bitprec);
2340 6412 : eno = get_eno(R,k,t,vi,v, vx, bitprec, 0);
2341 6412 : if (!eno && !thetad)
2342 : { /* t = sqrt(2), vi = theta(1/t), v = theta(t) */
2343 28 : lfunthetaspec(linit, bitprec, &vi, &v);
2344 28 : t = sqrtr(utor(2, prec));
2345 28 : eno = get_eno(R,k,t,conj_i(v),vi, vx, bitprec, 0);
2346 : }
2347 6412 : av = avma;
2348 6454 : while (!eno)
2349 : {
2350 42 : t = addsr(1, shiftr(utor(pari_rand(), prec), -2-BITS_IN_LONG));
2351 : /* t in [1,1.25[ */
2352 0 : v = thetad? lfuntheta(thetad, t, 0, bitprec)
2353 42 : : conj_i(lfuntheta(linit, t, 0, bitprec));
2354 42 : vi = lfuntheta(linit, ginv(t), 0, bitprec);
2355 42 : eno = get_eno(R,k,t,v,vi, vx, bitprec, c++ == 5);
2356 42 : set_avma(av);
2357 : }
2358 6412 : delete_var(); return ropm1(eno,prec);
2359 : }
2360 :
2361 : /* Find root number and/or residues when L-function coefficients and
2362 : conductor are known. For the moment at most a single residue allowed. */
2363 : GEN
2364 7364 : lfunrootres(GEN data, long bitprec)
2365 : {
2366 7364 : pari_sp ltop = avma;
2367 : GEN k, w, r, R, a, b, e, v, v2, be, ldata, linit;
2368 : long prec;
2369 :
2370 7364 : ldata = lfunmisc_to_ldata_shallow(data);
2371 7364 : r = ldata_get_residue(ldata);
2372 7364 : k = ldata_get_k(ldata);
2373 7364 : w = ldata_get_rootno(ldata);
2374 7364 : if (r) r = normalize_simple_pole(r, k);
2375 7364 : if (!r || residues_known(r))
2376 : {
2377 6951 : if (isintzero(w)) w = lfunrootno(data, bitprec);
2378 6951 : if (!r)
2379 4690 : r = R = gen_0;
2380 : else
2381 2261 : R = lfunrtoR_eno(ldata, w, nbits2prec(bitprec));
2382 6951 : return gc_GEN(ltop, mkvec3(r, R, w));
2383 : }
2384 413 : linit = lfunthetacheckinit(data, dbltor(sqrt(0.5)), 0, bitprec);
2385 413 : prec = nbits2prec(bitprec);
2386 413 : if (lg(r) > 2) pari_err_IMPL("multiple poles in lfunrootres");
2387 : /* Now residue unknown, and r = [[be,0]]. */
2388 413 : be = gmael(r, 1, 1);
2389 413 : if (ldata_isreal(ldata) && gequalm1(w))
2390 0 : R = lfuntheta(linit, gen_1, 0, bitprec);
2391 : else
2392 : {
2393 413 : GEN p2k = gpow(gen_2,k,prec);
2394 413 : lfunthetaspec(linit, bitprec, &v2, &v);
2395 413 : if (gequal(gmulsg(2, be), k)) pari_err_IMPL("pole at k/2 in lfunrootres");
2396 413 : if (gequal(be, k))
2397 : {
2398 147 : a = conj_i(gsub(gmul(p2k, v), v2));
2399 147 : b = subiu(p2k, 1);
2400 147 : e = gmul(gsqrt(p2k, prec), gsub(v2, v));
2401 : }
2402 : else
2403 : {
2404 266 : GEN tk2 = gsqrt(p2k, prec);
2405 266 : GEN tbe = gpow(gen_2, be, prec);
2406 266 : GEN tkbe = gpow(gen_2, gdivgu(gsub(k, be), 2), prec);
2407 266 : a = conj_i(gsub(gmul(tbe, v), v2));
2408 266 : b = gsub(gdiv(tbe, tkbe), tkbe);
2409 266 : e = gsub(gmul(gdiv(tbe, tk2), v2), gmul(tk2, v));
2410 : }
2411 413 : if (isintzero(w))
2412 : { /* Now residue unknown, r = [[be,0]], and w unknown. */
2413 7 : GEN t0 = mkfrac(utoi(11),utoi(10));
2414 7 : GEN th1 = lfuntheta(linit, t0, 0, bitprec);
2415 7 : GEN th2 = lfuntheta(linit, ginv(t0), 0, bitprec);
2416 7 : GEN tbe = gpow(t0, gmulsg(2, be), prec);
2417 7 : GEN tkbe = gpow(t0, gsub(k, be), prec);
2418 7 : GEN tk2 = gpow(t0, k, prec);
2419 7 : GEN c = conj_i(gsub(gmul(tbe, th1), th2));
2420 7 : GEN d = gsub(gdiv(tbe, tkbe), tkbe);
2421 7 : GEN f = gsub(gmul(gdiv(tbe, tk2), th2), gmul(tk2, th1));
2422 7 : GEN D = gsub(gmul(a, d), gmul(b, c));
2423 7 : w = gdiv(gsub(gmul(d, e), gmul(b, f)), D);
2424 : }
2425 413 : w = ropm1(w, prec);
2426 413 : R = gdiv(gsub(e, gmul(a, w)), b);
2427 : }
2428 413 : r = normalize_simple_pole(Rtor(be, R, ldata, prec), be);
2429 413 : R = lfunrtoR_i(ldata, r, w, prec);
2430 413 : return gc_GEN(ltop, mkvec3(r, R, w));
2431 : }
2432 :
2433 : /*******************************************************************/
2434 : /* Zeros */
2435 : /*******************************************************************/
2436 : struct lhardyz_t {
2437 : long bitprec, prec;
2438 : GEN linit;
2439 : };
2440 :
2441 : static GEN
2442 18563 : lfunhardyzeros(void *E, GEN t)
2443 : {
2444 18563 : struct lhardyz_t *S = (struct lhardyz_t*)E;
2445 18563 : GEN z = gprec_wensure(lfunhardy(S->linit, t, S->bitprec), S->prec);
2446 18563 : return typ(z) == t_VEC ? RgV_prod(z): z;
2447 : }
2448 :
2449 : /* initialize for computation on critical line up to height h, zero
2450 : * of order <= m */
2451 : static GEN
2452 581 : lfuncenterinit(GEN lmisc, double h, long m, long bitprec)
2453 : {
2454 581 : GEN ldata = lfunmisc_to_ldata_shallow(lmisc);
2455 581 : if (m < 0)
2456 : { /* choose a sensible default */
2457 581 : m = 4;
2458 581 : if (is_linit(lmisc) && linit_get_type(lmisc) == t_LDESC_INIT)
2459 : {
2460 490 : GEN domain = lfun_get_domain(linit_get_tech(lmisc));
2461 490 : m = domain_get_der(domain);
2462 : }
2463 : }
2464 581 : if (is_dirichlet(ldata)) m = 0;
2465 581 : return lfuninit(lmisc, mkvec(dbltor(h)), m, bitprec);
2466 : }
2467 :
2468 : long
2469 574 : lfunorderzero(GEN lmisc, long m, long bitprec)
2470 : {
2471 574 : pari_sp ltop = avma;
2472 : GEN eno, ldata, linit, k2;
2473 : long G, c0, c, st;
2474 :
2475 574 : if (is_linit(lmisc) && linit_get_type(lmisc) == t_LDESC_PRODUCT)
2476 : {
2477 84 : GEN M = gmael(linit_get_tech(lmisc), 2,1);
2478 84 : long i, l = lg(M);
2479 280 : for (c=0, i=1; i < l; i++) c += lfunorderzero(gel(M,i), m, bitprec);
2480 84 : return c;
2481 : }
2482 490 : linit = lfuncenterinit(lmisc, 0, m, bitprec);
2483 490 : ldata = linit_get_ldata(linit);
2484 490 : eno = ldata_get_rootno(ldata);
2485 490 : k2 = gmul2n(ldata_get_k(ldata), -1);
2486 490 : G = -bitprec/2;
2487 490 : c0 = 0; st = 1;
2488 490 : if (typ(eno) == t_VEC)
2489 : {
2490 42 : long i, l = lg(eno), cnt = l-1, s = 0;
2491 42 : GEN v = zero_zv(l-1);
2492 42 : if (ldata_isreal(ldata)) st = 2;
2493 84 : for (c = c0; cnt; c += st)
2494 : {
2495 42 : GEN L = lfun0(linit, k2, c, bitprec);
2496 154 : for (i = 1; i < l; i++)
2497 : {
2498 112 : if (v[i]==0 && gexpo(gel(L,i)) > G)
2499 : {
2500 112 : v[i] = c; cnt--; s += c;
2501 : }
2502 : }
2503 : }
2504 42 : return gc_long(ltop,s);
2505 : }
2506 : else
2507 : {
2508 448 : if (ldata_isreal(ldata)) { st = 2; if (!gequal1(eno)) c0 = 1; }
2509 448 : for (c = c0;; c += st)
2510 476 : if (gexpo(lfun0(linit, k2, c, bitprec)) > G) return gc_long(ltop, c);
2511 : }
2512 : }
2513 :
2514 : /* assume T1 * T2 > 0, T1 <= T2 */
2515 : static void
2516 98 : lfunzeros_i(struct lhardyz_t *S, GEN *pw, long *ct, GEN T1, GEN T2, long d,
2517 : GEN cN, GEN pi2, GEN pi2div, long precinit, long prec)
2518 : {
2519 98 : GEN T = T1, w = *pw;
2520 98 : long W = lg(w)-1, s = gsigne(lfunhardyzeros(S, T1));
2521 : for(;;)
2522 427 : {
2523 525 : pari_sp av = avma;
2524 : GEN D, T0, z;
2525 525 : D = gcmp(T, pi2) < 0? cN
2526 525 : : gadd(cN, gmulsg(d, glog(gdiv(T, pi2), prec)));
2527 525 : D = gdiv(pi2div, gmulsg(d, D));
2528 : for(;;)
2529 13482 : {
2530 : long s0;
2531 14007 : T0 = T; T = gadd(T, D);
2532 14007 : if (gcmp(T, T2) >= 0) T = T2;
2533 14007 : s0 = gsigne(lfunhardyzeros(S, T));
2534 14007 : if (s0 != s) { s = s0; break; }
2535 13580 : if (T == T2) { setlg(w, *ct); *pw = w; return; }
2536 : }
2537 427 : z = zbrent(S, lfunhardyzeros, T0, T, prec); /* T <= T2 */
2538 427 : (void)gc_all(av, 2, &T, &z);
2539 427 : if (*ct > W) { W *= 2; w = vec_lengthen(w, W); }
2540 427 : if (typ(z) == t_REAL) z = rtor(z, precinit);
2541 427 : gel(w, (*ct)++) = z;
2542 : }
2543 : setlg(w, *ct); *pw = w;
2544 : }
2545 : GEN
2546 98 : lfunzeros(GEN ldata, GEN lim, long divz, long bitprec)
2547 : {
2548 98 : pari_sp ltop = avma;
2549 : GEN linit, pi2, pi2div, cN, w, T, h1, h2;
2550 98 : long i, d, NEWD, c, ct, s1, s2, prec, prec0 = nbits2prec(bitprec);
2551 : double maxt;
2552 : struct lhardyz_t S;
2553 :
2554 98 : if (is_linit(ldata) && linit_get_type(ldata) == t_LDESC_PRODUCT)
2555 : {
2556 0 : GEN M = gmael(linit_get_tech(ldata), 2,1);
2557 0 : long l = lg(M);
2558 0 : w = cgetg(l, t_VEC);
2559 0 : for (i = 1; i < l; i++) gel(w,i) = lfunzeros(gel(M,i), lim, divz, bitprec);
2560 0 : return gc_upto(ltop, vecsort0(shallowconcat1(w), NULL, 0));
2561 : }
2562 98 : if (typ(lim) == t_VEC)
2563 : {
2564 : double H1, H2;
2565 63 : if (lg(lim) != 3 || gcmp(gel(lim, 1), gel(lim, 2)) >= 0)
2566 7 : pari_err_TYPE("lfunzeros",lim);
2567 56 : h1 = gel(lim, 1); H1 = gtodouble(h1);
2568 56 : h2 = gel(lim, 2); H2 = gtodouble(h2);
2569 56 : maxt = maxdd(fabs(H1), fabs(H2));
2570 56 : if (H1 * H2 > 0)
2571 : {
2572 35 : GEN LDATA = lfunmisc_to_ldata_shallow(ldata);
2573 35 : double m = mindd(fabs(H1), fabs(H2));
2574 35 : if (is_dirichlet(LDATA) && m > lfuninit_cutoff(LDATA)) maxt = 0;
2575 : }
2576 : }
2577 : else
2578 : {
2579 35 : if (gcmp(lim, gen_0) <= 0) pari_err_TYPE("lfunzeros",lim);
2580 35 : h1 = gen_0;
2581 35 : h2 = lim;
2582 35 : maxt = gtodouble(h2);
2583 : }
2584 91 : S.linit = linit = lfuncenterinit(ldata, maxt, -1, bitprec);
2585 91 : S.bitprec = bitprec;
2586 91 : S.prec = prec0;
2587 91 : ldata = linit_get_ldata(linit);
2588 91 : d = ldata_get_degree(ldata);
2589 :
2590 91 : NEWD = minss((long) ceil(bitprec + (M_PI/(4*M_LN2)) * d * maxt),
2591 : lfun_get_bitprec(linit_get_tech(linit)));
2592 91 : prec = nbits2prec(NEWD);
2593 91 : cN = gdiv(ldata_get_conductor(ldata), gpowgs(Pi2n(-1, prec), d));
2594 91 : cN = gexpo(cN) >= 0? gaddsg(d, gmulsg(2, glog(cN, prec))): utoi(d);
2595 91 : pi2 = Pi2n(1, prec);
2596 91 : pi2div = gdivgu(pi2, labs(divz));
2597 91 : s1 = gsigne(h1);
2598 91 : s2 = gsigne(h2);
2599 91 : w = cgetg(100+1, t_VEC); c = 1; ct = 0; T = NULL;
2600 91 : if (s1 <= 0 && s2 >= 0)
2601 : {
2602 56 : GEN r = ldata_get_residue(ldata);
2603 56 : if (!r || gequal0(r))
2604 : {
2605 35 : ct = lfunorderzero(linit, -1, bitprec);
2606 35 : if (ct) T = real2n(-prec / (2*ct), prec);
2607 : }
2608 : }
2609 91 : if (s1 <= 0)
2610 : {
2611 63 : if (s1 < 0)
2612 21 : lfunzeros_i(&S, &w, &c, h1, T? negr(T): h2,
2613 : d, cN, pi2, pi2div, prec0, prec);
2614 63 : if (ct)
2615 : {
2616 21 : long n = lg(w)-1;
2617 21 : if (c + ct >= n) w = vec_lengthen(w, n + ct);
2618 84 : for (i = 1; i <= ct; i++) gel(w,c++) = gen_0;
2619 : }
2620 : }
2621 91 : if (s2 > 0 && (T || s1 >= 0))
2622 77 : lfunzeros_i(&S, &w, &c, T? T: h1, h2, d, cN, pi2, pi2div, prec0, prec);
2623 91 : return gc_GEN(ltop, w);
2624 : }
2625 :
2626 : /*******************************************************************/
2627 : /* Guess conductor */
2628 : /*******************************************************************/
2629 : struct huntcond_t {
2630 : GEN k;
2631 : GEN theta, thetad;
2632 : GEN *pM, *psqrtM, *pMd, *psqrtMd;
2633 : };
2634 :
2635 : static void
2636 11879 : condset(struct huntcond_t *S, GEN M, long prec)
2637 : {
2638 11879 : *(S->pM) = M;
2639 11879 : *(S->psqrtM) = gsqrt(ginv(M), prec);
2640 11879 : if (S->thetad != S->theta)
2641 : {
2642 0 : *(S->pMd) = *(S->pM);
2643 0 : *(S->psqrtMd) = *(S->psqrtM);
2644 : }
2645 11879 : }
2646 :
2647 : /* M should eventually converge to N, the conductor. L has no pole. */
2648 : static GEN
2649 6895 : wrap1(void *E, GEN M)
2650 : {
2651 6895 : struct huntcond_t *S = (struct huntcond_t*)E;
2652 : GEN thetainit, tk, p1, p1inv;
2653 6895 : GEN t = mkfrac(stoi(11), stoi(10));
2654 : long prec, bitprec;
2655 :
2656 6895 : thetainit = linit_get_tech(S->theta);
2657 6895 : bitprec = theta_get_bitprec(thetainit);
2658 6895 : prec = nbits2prec(bitprec);
2659 6895 : condset(S, M, prec);
2660 6895 : tk = gpow(t, S->k, prec);
2661 6895 : p1 = lfuntheta(S->thetad, t, 0, bitprec);
2662 6895 : p1inv = lfuntheta(S->theta, ginv(t), 0, bitprec);
2663 6895 : return glog(gabs(gmul(tk, gdiv(p1, p1inv)), prec), prec);
2664 : }
2665 :
2666 : /* M should eventually converge to N, the conductor. L has a pole. */
2667 : static GEN
2668 4942 : wrap2(void *E, GEN M)
2669 : {
2670 4942 : struct huntcond_t *S = (struct huntcond_t*)E;
2671 : GEN t1k, t2k, p1, p1inv, p2, p2inv, thetainit, R;
2672 4942 : GEN t1 = mkfrac(stoi(11), stoi(10)), t2 = mkfrac(stoi(13), stoi(11));
2673 : GEN t1be, t2be, t1bemk, t2bemk, t1kmbe, t2kmbe;
2674 : GEN F11, F12, F21, F22, P1, P2, res;
2675 : long prec, bitprec;
2676 4942 : GEN k = S->k;
2677 :
2678 4942 : thetainit = linit_get_tech(S->theta);
2679 4942 : bitprec = theta_get_bitprec(thetainit);
2680 4942 : prec = nbits2prec(bitprec);
2681 4942 : condset(S, M, prec);
2682 :
2683 4942 : p1 = lfuntheta(S->thetad, t1, 0, bitprec);
2684 4942 : p2 = lfuntheta(S->thetad, t2, 0, bitprec);
2685 4942 : p1inv = lfuntheta(S->theta, ginv(t1), 0, bitprec);
2686 4942 : p2inv = lfuntheta(S->theta, ginv(t2), 0, bitprec);
2687 4942 : t1k = gpow(t1, k, prec);
2688 4942 : t2k = gpow(t2, k, prec);
2689 4942 : R = theta_get_R(thetainit);
2690 4942 : if (typ(R) == t_VEC)
2691 : {
2692 0 : GEN be = gmael(R, 1, 1);
2693 0 : t1be = gpow(t1, be, prec); t1bemk = gdiv(gsqr(t1be), t1k);
2694 0 : t2be = gpow(t2, be, prec); t2bemk = gdiv(gsqr(t2be), t2k);
2695 0 : t1kmbe = gdiv(t1k, t1be);
2696 0 : t2kmbe = gdiv(t2k, t2be);
2697 : }
2698 : else
2699 : { /* be = k */
2700 4942 : t1be = t1k; t1bemk = t1k; t1kmbe = gen_1;
2701 4942 : t2be = t2k; t2bemk = t2k; t2kmbe = gen_1;
2702 : }
2703 4942 : F11 = conj_i(gsub(gmul(gsqr(t1be), p1), p1inv));
2704 4942 : F12 = conj_i(gsub(gmul(gsqr(t2be), p2), p2inv));
2705 4942 : F21 = gsub(gmul(t1k, p1), gmul(t1bemk, p1inv));
2706 4942 : F22 = gsub(gmul(t2k, p2), gmul(t2bemk, p2inv));
2707 4942 : P1 = gsub(gmul(t1bemk, t1be), t1kmbe);
2708 4942 : P2 = gsub(gmul(t2bemk, t2be), t2kmbe);
2709 4942 : res = gdiv(gsub(gmul(P2,F21), gmul(P1,F22)),
2710 : gsub(gmul(P2,F11), gmul(P1,F12)));
2711 4942 : return glog(gabs(res, prec), prec);
2712 : }
2713 :
2714 : /* If flag = 0 (default) return all conductors found as integers. If
2715 : flag = 1, return the approximations, not the integers. If flag = 2,
2716 : return all, even nonintegers. */
2717 :
2718 : static GEN
2719 84 : checkconductor(GEN v, long bit, long flag)
2720 : {
2721 : GEN w;
2722 84 : long e, j, k, l = lg(v);
2723 84 : if (flag == 2) return v;
2724 84 : w = cgetg(l, t_VEC);
2725 322 : for (j = k = 1; j < l; j++)
2726 : {
2727 238 : GEN N = grndtoi(gel(v,j), &e);
2728 238 : if (e < -bit) gel(w,k++) = flag ? gel(v,j): N;
2729 : }
2730 84 : if (k == 2) return gel(w,1);
2731 7 : setlg(w,k); return w;
2732 : }
2733 :
2734 : static GEN
2735 98 : parse_maxcond(GEN maxN)
2736 : {
2737 : GEN M;
2738 98 : if (!maxN)
2739 49 : M = utoipos(10000);
2740 49 : else if (typ(maxN) == t_VEC)
2741 : {
2742 14 : if (!RgV_is_ZV(maxN)) pari_err_TYPE("lfunconductor",maxN);
2743 14 : return ZV_sort_shallow(maxN);
2744 : }
2745 : else
2746 35 : M = maxN;
2747 84 : return (typ(M) == t_INT)? addiu(M, 1): gceil(M);
2748 : }
2749 :
2750 : GEN
2751 98 : lfunconductor(GEN data, GEN maxcond, long flag, long bitprec)
2752 : {
2753 : struct huntcond_t S;
2754 98 : pari_sp av = avma;
2755 98 : GEN ldata = lfunmisc_to_ldata_shallow(data);
2756 98 : GEN ld, r, v, theta, thetad, M, tdom, t0 = NULL, t0i = NULL;
2757 : GEN (*eval)(void *, GEN);
2758 : long prec;
2759 98 : M = parse_maxcond(maxcond);
2760 98 : r = ldata_get_residue(ldata);
2761 98 : if (typ(M) == t_VEC) /* select in list */
2762 : {
2763 14 : if (lg(M) == 1) retgc_const(av, cgetg(1, t_VEC));
2764 7 : eval = NULL; tdom = dbltor(0.7);
2765 : }
2766 84 : else if (!r) { eval = wrap1; tdom = uutoQ(10,11); }
2767 : else
2768 : {
2769 21 : if (typ(r) == t_VEC && lg(r) > 2)
2770 0 : pari_err_IMPL("multiple poles in lfunconductor");
2771 21 : eval = wrap2; tdom = uutoQ(11,13);
2772 : }
2773 91 : if (eval) bitprec += bitprec/2;
2774 91 : prec = nbits2prec(bitprec);
2775 91 : ld = shallowcopy(ldata);
2776 91 : gel(ld, 5) = eval? M: veclast(M);
2777 91 : theta = lfunthetainit_i(ld, tdom, 0, bitprec);
2778 91 : thetad = theta_dual(theta, ldata_get_dual(ldata));
2779 91 : gel(theta,3) = shallowcopy(linit_get_tech(theta));
2780 91 : S.k = ldata_get_k(ldata);
2781 91 : S.theta = theta;
2782 91 : S.thetad = thetad? thetad: theta;
2783 91 : S.pM = &gel(linit_get_ldata(theta),5);
2784 91 : S.psqrtM = &gel(linit_get_tech(theta),7);
2785 91 : if (thetad)
2786 : {
2787 0 : S.pMd = &gel(linit_get_ldata(thetad),5);
2788 0 : S.psqrtMd = &gel(linit_get_tech(thetad),7);
2789 : }
2790 91 : if (!eval)
2791 : {
2792 7 : long i, besti = 0, beste = -10, l = lg(M);
2793 7 : t0 = uutoQ(11,10); t0i = uutoQ(10,11);
2794 49 : for (i = 1; i < l; i++)
2795 : {
2796 42 : pari_sp av2 = avma;
2797 : long e;
2798 42 : condset(&S, gel(M,i), prec);
2799 42 : e = lfuncheckfeq_i(theta, thetad, t0, t0i, bitprec);
2800 42 : set_avma(av2);
2801 42 : if (e < beste) { beste = e; besti = i; }
2802 35 : else if (e == beste) beste = besti = 0; /* tie: forget */
2803 : }
2804 7 : if (!besti) retgc_const(av, cgetg(1, t_VEC));
2805 7 : return gc_GEN(av, mkvec2(gel(M,besti), stoi(beste)));
2806 : }
2807 84 : v = solvestep((void*)&S, eval, ghalf, M, gen_2, 14, prec);
2808 84 : return gc_GEN(av, checkconductor(v, bitprec/2, flag));
2809 : }
2810 :
2811 : /* assume chi primitive */
2812 : static GEN
2813 2863 : znchargauss_i(GEN G, GEN chi, long bitprec)
2814 : {
2815 2863 : GEN z, q, F = znstar_get_N(G);
2816 : long prec;
2817 :
2818 2863 : if (equali1(F)) return gen_1;
2819 2863 : prec = nbits2prec(bitprec);
2820 2863 : q = sqrtr_abs(itor(F, prec));
2821 2863 : z = lfuntheta(mkvec2(G,chi), gen_1, 0, bitprec);
2822 2863 : if (gexpo(z) < 10 - bitprec)
2823 : {
2824 28 : if (equaliu(F,300))
2825 : {
2826 14 : GEN z = rootsof1u_cx(25, prec);
2827 14 : GEN n = znconreyexp(G, chi);
2828 14 : if (equaliu(n, 131)) return gmul(q, gpowgs(z,14));
2829 7 : if (equaliu(n, 71)) return gmul(q, gpowgs(z,11));
2830 : }
2831 14 : if (equaliu(F,600))
2832 : {
2833 14 : GEN z = rootsof1u_cx(25, prec);
2834 14 : GEN n = znconreyexp(G, chi);
2835 14 : if (equaliu(n, 491)) return gmul(q, gpowgs(z,7));
2836 7 : if (equaliu(n, 11)) return gmul(q, gpowgs(z,18));
2837 : }
2838 0 : pari_err_BUG("znchargauss [ Theta(chi,1) = 0 ]");
2839 : }
2840 2835 : z = gmul(gdiv(z, conj_i(z)), q);
2841 2835 : if (zncharisodd(G,chi)) z = mulcxI(z);
2842 2835 : return z;
2843 : }
2844 : static GEN
2845 2863 : Z_radical(GEN N, long *om)
2846 : {
2847 2863 : GEN P = gel(Z_factor(N), 1);
2848 2863 : *om = lg(P)-1; return ZV_prod(P);
2849 : }
2850 : GEN
2851 5516 : znchargauss(GEN G, GEN chi, GEN a, long bitprec)
2852 : {
2853 : GEN v, T, N, F, b0, b1, b2, bF, a1, aF, A, r, GF, tau, B, faB, u, S;
2854 5516 : long omb0, prec = nbits2prec(bitprec);
2855 5516 : pari_sp av = avma;
2856 :
2857 5516 : if (typ(chi) != t_COL) chi = znconreylog(G,chi);
2858 5516 : T = znchartoprimitive(G, chi);
2859 5516 : GF = gel(T,1);
2860 5516 : chi = gel(T,2); /* now primitive */
2861 5516 : N = znstar_get_N(G);
2862 5516 : F = znstar_get_N(GF);
2863 5516 : if (equalii(N,F)) b1 = bF = gen_1;
2864 : else
2865 : {
2866 245 : v = Z_ppio(diviiexact(N,F), F);
2867 245 : bF = gel(v,2); /* (N/F, F^oo) */
2868 245 : b1 = gel(v,3); /* cofactor */
2869 : }
2870 5516 : if (!a) a = a1 = aF = gen_1;
2871 : else
2872 : {
2873 5467 : if (typ(a) != t_INT) pari_err_TYPE("znchargauss",a);
2874 5467 : a = modii(a, N);
2875 5467 : if (!signe(a)) { set_avma(av); return is_pm1(F)? eulerphi(N): gen_0; }
2876 3031 : v = Z_ppio(a, F);
2877 3031 : aF = gel(v,2);
2878 3031 : a1 = gel(v,3);
2879 : }
2880 3080 : if (!equalii(aF, bF)) { set_avma(av); return gen_0; }
2881 2863 : b0 = Z_radical(b1, &omb0);
2882 2863 : b2 = diviiexact(b1, b0);
2883 2863 : A = dvmdii(a1, b2, &r);
2884 2863 : if (r != gen_0) { set_avma(av); return gen_0; }
2885 2863 : B = gcdii(A,b0); faB = Z_factor(B); /* squarefree */
2886 2863 : S = eulerphi(mkvec2(B,faB));
2887 2863 : if (odd(omb0 + lg(gel(faB,1))-1)) S = negi(S); /* moebius(b0/B) * phi(B) */
2888 2863 : S = mulii(S, mulii(aF,b2));
2889 2863 : tau = znchargauss_i(GF, chi, bitprec);
2890 2863 : u = Fp_div(b0, A, F);
2891 2863 : if (!equali1(u))
2892 : {
2893 7 : GEN ord = zncharorder(GF, chi), z = rootsof1_cx(ord, prec);
2894 7 : tau = gmul(tau, znchareval(GF, chi, u, mkvec2(z,ord)));
2895 : }
2896 2863 : return gc_upto(av, gmul(tau, S));
2897 : }
|
