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