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 : #include "pari.h"
16 : #include "paripriv.h"
17 :
18 : #define DEBUGLEVEL DEBUGLEVEL_gammamellininv
19 :
20 : /*******************************************************************/
21 : /* Computation of inverse Mellin */
22 : /* transforms of gamma products. */
23 : /*******************************************************************/
24 : /* Handle complex Vga whose sum is real */
25 : static GEN
26 24514 : sumVga(GEN Vga) { return real_i(vecsum(Vga)); }
27 :
28 : /* ac != 0 */
29 : static double
30 568524 : lemma526_i(double ac, double c, double t, double B)
31 : {
32 568524 : double D = -B/ac; /* sgn(t) = sgn(a) = - sgn(D) */
33 568524 : if (D <= 0)
34 : {
35 493944 : if (D > -100)
36 : {
37 67298 : D = -exp(D) / t;
38 67298 : if (D < - 1/M_E) return 0;
39 66969 : D = dbllambertW_1(D);
40 : }
41 : else
42 : { /* avoid underflow, use asymptotic expansion */
43 426646 : double U = D - log(t);
44 426646 : D = U - log(-U);
45 : }
46 493615 : return pow(maxdd(t, -t * D), c);
47 : }
48 : else
49 : {
50 74580 : if (D < 100)
51 9856 : D = dbllambertW0(-exp(D) / t);
52 : else
53 : { /* avoid overflow, use asymptotic expansion */
54 64724 : double U = D - log(-t);
55 64724 : D = U - log(U);
56 : }
57 74580 : return pow(-t * D, c);
58 : }
59 : }
60 : /* b > 0, c > 0; solve x^a exp(-b x^(1/c)) < e^(-B) for x >= 0 */
61 : double
62 14 : dbllemma526(double a, double b, double c, double B)
63 : {
64 : double ac;
65 14 : if (!a) return B <= 0? 0: pow(B/b, c);
66 14 : ac = a*c; if (B < 0) B = 1e-9;
67 14 : return lemma526_i(ac, c, ac/b, B);
68 : }
69 : /* Same, special case b/c = 2Pi, the only one needed: for c = d/2 */
70 : double
71 2503046 : dblcoro526(double a, double c, double B)
72 : {
73 2503046 : if (!a) return B <= 0? 0: pow(B/(2*M_PI*c), c);
74 568510 : if (B < 0) B = 1e-9;
75 568510 : return lemma526_i(a*c, c, a/(2*M_PI), B);
76 : }
77 :
78 : static const double MELLININV_CUTOFF = 121.; /* C*C */
79 :
80 : /* x real */
81 : static GEN
82 9506 : RMOD2(GEN x) { return gsub(x, gmul2n(gdiventgs(x,2), 1)); }
83 : /* x real or complex, return canonical representative for x mod 2Z */
84 : static GEN
85 9506 : MOD2(GEN x)
86 9506 : { return typ(x) == t_COMPLEX? mkcomplex(RMOD2(gel(x,1)), gel(x,2)): RMOD2(x); }
87 : static GEN
88 2912 : RgV_MOD2(GEN x)
89 12418 : { pari_APPLY_same(MOD2(gel(x,i))); }
90 :
91 : /* classes of poles of the gamma factor mod 2Z, sorted by increasing
92 : * Re(s) mod 2 (in [0,2[).*/
93 : static GEN
94 2912 : gammapoles(GEN Vga, long *pdV, long bit)
95 : {
96 2912 : long i, m, emax, l = lg(Vga);
97 2912 : GEN P, B = RgV_MOD2(Vga), V = cgetg(l, t_VEC);
98 2912 : P = gen_indexsort(B, (void*)lexcmp, cmp_nodata);
99 8036 : for (i = m = 1; i < l;)
100 : {
101 5124 : GEN u = gel(B, P[i]);
102 : long k;
103 9506 : for(k = i+1; k < l; k++)
104 : {
105 6594 : GEN v = gsub(u, gel(B, P[k]));
106 6594 : if (!gequal0(v) && (!isinexactreal(v) || gexpo(v) > -bit)) break;
107 : }
108 5124 : gel(V, m++) = vecslice(P,i,k-1);
109 5124 : i = k;
110 : }
111 2912 : setlg(V, m); emax = 0;
112 8036 : for (i = 1; i < m; i++)
113 : {
114 5124 : long j, e = 0, li = lg(gel(V,i))-1;
115 5124 : GEN b = gel(B, gel(V,i)[1]);
116 15610 : for (j = 1; j < m; j++)
117 10486 : if (j != i) e -= gexpo(gsub(gel(B, gel(V,j)[1]), b));
118 5124 : emax = maxss(emax, e * li);
119 : }
120 8036 : for (i = 1; i < m; i++) gel(V,i) = vecpermute(Vga, gel(V,i));
121 2912 : *pdV = emax; return V;
122 : }
123 :
124 : static GEN
125 476644 : sercoeff(GEN x, long n, long prec)
126 : {
127 476644 : long N = n - valser(x);
128 476644 : return (N < 0)? gen_0: gprec_wtrunc(gel(x, N+2), prec);
129 : }
130 :
131 : /* prod_i Gamma(s/2 + (m+LA[i])/2), set t *= prod_i (s/2 + (m+LA[i])/2) */
132 : static GEN
133 10486 : get_gamma(GEN *pt, GEN LA, GEN m, int round, long precdl, long prec)
134 : {
135 10486 : long i, l = lg(LA);
136 10486 : GEN pr = NULL, t = *pt;
137 27979 : for (i = 1; i < l; i++)
138 : {
139 17493 : GEN u, g, a = gmul2n(gadd(m, gel(LA,i)), -1);
140 17493 : if (round) a = ground(a);
141 17493 : u = deg1pol_shallow(ghalf, a, 0);
142 17493 : g = ggamma(RgX_to_ser(u, precdl), prec);
143 17493 : pr = pr? gmul(pr, g): g;
144 17493 : t = t? gmul(t, u): u;
145 : }
146 10486 : *pt = t; return pr;
147 : }
148 : /* generalized power series expansion of inverse Mellin around x = 0;
149 : * m-th derivative */
150 : static GEN
151 2912 : Kderivsmallinit(GEN ldata, GEN Vga, long m, long bit)
152 : {
153 2912 : const double C2 = MELLININV_CUTOFF;
154 2912 : long prec2, N, j, l, dLA, limn, d = lg(Vga)-1;
155 : GEN piA, LA, L, M, mat;
156 :
157 2912 : LA = gammapoles(Vga, &dLA, bit); N = lg(LA)-1;
158 2912 : prec2 = nbits2prec(dLA + bit * (1 + M_PI*d/C2));
159 : #if BITS_IN_LONG == 32
160 416 : prec2 += odd(prec2lg(prec2)) ? EXTRAPRECWORD: 0;
161 : #endif
162 2912 : if (ldata) Vga = ldata_get_gammavec(ldata_newprec(ldata, prec2));
163 2912 : L = cgetg(N+1, t_VECSMALL);
164 2912 : M = cgetg(N+1, t_VEC);
165 2912 : mat = cgetg(N+1, t_VEC);
166 2912 : limn = ceil(2*M_LN2*bit / (d * dbllambertW0(C2/(M_PI*M_E*d))));
167 2912 : l = limn + 2;
168 8036 : for (j = 1; j <= N; j++)
169 : {
170 5124 : GEN S = gel(LA,j);
171 5124 : GEN C, c, mj, G = NULL, t = NULL, tj = NULL;
172 5124 : long i, k, n, jj, lj = L[j] = lg(S)-1, precdl = lj+3;
173 :
174 5124 : gel(M,j) = mj = gsubsg(2, gel(S, vecindexmin(real_i(S))));
175 15610 : for (jj = 1; jj <= N; jj++)
176 : {
177 : GEN g;
178 10486 : if (jj == j) /* poles come from this class only */
179 5124 : g = get_gamma(&tj, gel(LA,jj), mj, 1, precdl, prec2);
180 : else
181 5362 : g = get_gamma(&t, gel(LA,jj), mj, 0, precdl, prec2);
182 10486 : G = G? gmul(G, g): g;
183 : }
184 5124 : c = cgetg(limn+2,t_COL); gel(c,1) = G;
185 273042 : for (n=1; n <= limn; n++)
186 : {
187 267918 : GEN A = utoineg(2*n), T = RgX_translate(tj, A);
188 : /* T = exact polynomial, may vanish at 0 (=> pole in c[n+1]) */
189 267918 : if (t) T = RgX_mul(T, RgX_translate(t, A)); /* no pole here */
190 267918 : gel(c,n+1) = gdiv(gel(c,n), T);
191 : }
192 5124 : gel(mat, j) = C = cgetg(lj+1, t_COL);
193 14630 : for (k = 1; k <= lj; k++)
194 : {
195 9506 : GEN L = cgetg(l, t_POL);
196 486150 : for (n = 2; n < l; n++) gel(L,n) = sercoeff(gel(c,n), -k, prec2);
197 9506 : L[1] = evalsigne(1)|evalvarn(0); gel(C,k) = L;
198 : }
199 : /* C[k] = \sum_n c_{j,k} t^n =: C_k(t) in Dokchitser's Algo 3.3
200 : * m-th derivative of t^(-M+2) sum_k (-ln t)^k/k! C_k(t^2) */
201 5124 : if (m)
202 : {
203 119 : mj = gsubgs(mj, 2);
204 238 : for (i = 1; i <= m; i++, mj = gaddgs(mj,1))
205 322 : for (k = 1; k <= lj; k++)
206 : {
207 203 : GEN c = gel(C,k), d = RgX_shift_shallow(gmul2n(RgX_deriv(c),1), 1);
208 203 : c = RgX_Rg_mul(c, mj);
209 203 : if (k < lj) c = RgX_add(c, gel(C,k+1));
210 203 : gel(C,k) = RgX_sub(d, c);
211 : }
212 119 : gel(M,j) = gaddgs(mj,2);
213 : }
214 14630 : for (k = 1; k <= lj; k++)
215 : {
216 9506 : GEN c = gel(C,k);
217 9506 : if (k > 2) c = RgX_Rg_div(c, mpfact(k-1));
218 9506 : gel(C,k) = RgX_to_RgC(c, lgpol(c));
219 : }
220 : }
221 : /* Algo 3.3: * \phi^(m)(t) = sum_j t^m_j sum_k (-ln t)^k mat[j,k](t^2) */
222 2912 : piA = gsubsg(m*d, sumVga(Vga));
223 2912 : if (!gequal0(piA)) piA = powPis(gmul2n(piA,-1), prec2);
224 2912 : return mkvec5(L, RgV_neg(M), mat, mkvecsmall(prec2), piA);
225 : }
226 :
227 : /* Evaluate a vector considered as a polynomial using Horner. */
228 : static GEN
229 262284 : evalvec(GEN vec, long N, GEN u)
230 : {
231 262284 : GEN S = gen_0;
232 : long n;
233 262284 : N = minss(N, lg(vec)-1);
234 9111869 : for (n = N; n >= 1; n--) S = gmul(u, gadd(gel(vec,n), S));
235 262287 : return S;
236 : }
237 :
238 : /* gammamellininvinit accessors */
239 : static double
240 4522976 : get_tmax(long bitprec)
241 4522976 : { return (M_LN2 / MELLININV_CUTOFF) * bitprec ; }
242 : static GEN
243 5726 : GMi_get_Vga(GEN K) { return gel(K,2); }
244 : static long
245 11321136 : GMi_get_degree(GEN K) { return lg(gel(K,2))-1; }
246 : static long
247 5620403 : GMi_get_m(GEN K) { return itos( gel(K,3) ); }
248 : static GEN /* [lj,mj,mat,prec2], Kderivsmall only */
249 5771223 : GMi_get_VS(GEN K) { return gel(K,4); }
250 : /* K[5] = [Ms,cd,A2], Kderivlarge only */
251 : static long/*Kderivlarge*/
252 5620399 : GMi_get_status(GEN K) { return itos(gmael3(K,5,1,2)); }
253 : static GEN/*Kderivlarge*/
254 5620454 : GMi_get_M(GEN K) { return gmael3(K,5,1,1); }
255 : static GEN/*Kderivlarge*/
256 5620426 : GMi_get_cd(GEN K) { return gmael(K,5,2); }
257 : static GEN/*Kderivlarge*/
258 11239730 : GMi_get_A2(GEN K) { return gmael(K,5,3); }
259 :
260 : static double
261 5695951 : GMi_get_tmax(GEN K, long bitprec)
262 5695951 : { return (typ(GMi_get_VS(K)) == t_INT)? -1.0 : get_tmax(bitprec); }
263 :
264 : /* Compute m-th derivative of inverse Mellin at x by generalized power series
265 : * around x = 0; x2d = x^(2/d), x is possibly NULL (don't bother about
266 : * complex branches). Assume |x|^(2/d) <= tmax = M_LN2*bitprec/MELLININV_CUTOFF*/
267 : static GEN
268 75377 : Kderivsmall(GEN K, GEN x, GEN x2d, long bitprec)
269 : {
270 75377 : GEN VS = GMi_get_VS(K), L = gel(VS,1), M = gel(VS,2), mat = gel(VS,3);
271 75377 : GEN d2, Lx, x2, S, pi, piA = gel(VS,5);
272 75377 : long j, k, prec = gel(VS,4)[1], d = GMi_get_degree(K), limn, N = lg(L)-1;
273 75377 : double xd, Wd, Ed = M_LN2*bitprec / d;
274 :
275 75377 : xd = maxdd(M_PI*dblmodulus(x2d), 1E-13); /* pi |x|^2/d unless x tiny */
276 75377 : if (xd > Ed) pari_err_BUG("Kderivsmall (x2d too large)");
277 : /* Lemma 5.2.6 (2), a = 1 + log(Pi x^(2/d)) = log(e / xd),
278 : * B = log(2)*bitprec / d = Ed */
279 75377 : Wd = dbllambertW0( Ed / (M_E*xd) ); /* solution of w exp(w) = B exp(-a)*/
280 75377 : limn = (long) ceil(2*Ed/Wd);
281 75377 : pi = mppi(prec);
282 75377 : d2 = gdivsg(d,gen_2);
283 75377 : if (x)
284 1122 : x = gmul(gtofp(x,prec), gpow(pi,d2,prec));
285 : else
286 74255 : x = gpow(gmul(gtofp(x2d,prec),pi), d2, prec);
287 : /* at this stage, x has been replaced by pi^(d/2) x */
288 75376 : x2 = gsqr(x);
289 75377 : Lx = gpowers(gneg(glog(x,prec)), vecsmall_max(L));
290 75375 : S = gen_0;
291 205594 : for (j = 1; j <= N; ++j)
292 : {
293 130218 : long lj = L[j];
294 130218 : GEN s = gen_0;
295 392502 : for (k = 1; k <= lj; k++)
296 262284 : s = gadd(s, gmul(gel(Lx,k), evalvec(gmael(mat,j,k), limn, x2)));
297 130218 : S = gadd(S, gmul(gpow(x, gel(M,j), prec), s));
298 : }
299 75376 : if (!gequal0(piA)) S = gmul(S, piA);
300 75376 : return S;
301 : }
302 :
303 : /* In Klarge, we conpute K(t) as (asymptotic) * F(z), where F ~ 1 is given by
304 : * a continued fraction and z = Pi t^(2/d). If we take 2n terms in F (n terms
305 : * in Euler form), F_n(z) - F(z) is experimentally in exp(- C sqrt(n*z))
306 : * where C ~ 8 for d > 2 [HEURISTIC] and C = 4 (theorem) for d = 1 or d = 2
307 : * and vga = [0,1]. For e^(-E) absolute error, we want
308 : * exp(-C sqrt(nz)) < e^-(E+a), where a ~ ln(asymptotic)
309 : * i.e. 2n > (E+a)^2 / t^(2/d) * 2/(C^2 Pi); C^2*Pi/2 ~ 100.5 ~ 101
310 : *
311 : * In fact, this model becomes wrong for z large: we use instead
312 : *
313 : * exp(- sqrt(D * nz/log(z+1))) < e^-(E+a),
314 : * i.e. 2n > (E+a)^2 * log(1 + Pi t^(2/d))/ t^(2/d) * 2/(D Pi); */
315 : static double
316 5641780 : get_D(long d) { return d <= 2 ? 157. : 180.; }
317 : /* if (abs), absolute error rather than relative */
318 : static void
319 5620398 : Kderivlarge_optim(GEN K, int abs, GEN t2d, double cd, long *pbitprec, long *pnlim)
320 : {
321 5620398 : GEN A2 = GMi_get_A2(K);
322 5620360 : long bitprec = *pbitprec, d = GMi_get_degree(K);
323 5620428 : const double D = get_D(d), td = dblmodulus(t2d);
324 5620393 : double a, rtd, E = M_LN2*bitprec;
325 :
326 : /* t = 0 can happen with finite continued fraction or easyvga */
327 5620393 : if (!td) { *pnlim = 0; return; }
328 5620386 : rtd = (typ(t2d) == t_COMPLEX)? gtodouble(gel(t2d,1)): td;
329 : /* A2/2 = A, log(td) = (2/d)*log t */
330 5620386 : a = d*gtodouble(A2)*log2(td)/2 - (M_PI/M_LN2)*d*rtd + log2(cd);/*log2 K(t)~a*/
331 : /* if bitprec <= 0, caller should return K(t) ~ 0 */
332 5619936 : bitprec += 64;
333 5619936 : if (abs)
334 : {
335 5615247 : bitprec += ceil(a);
336 5615247 : if (a <= -65) E = M_LN2*bitprec; /* guarantees E <= initial E */
337 : }
338 5619936 : *pbitprec = bitprec;
339 5619936 : *pnlim = ceil(E*E * log2(1+M_PI*td) / (D*td));
340 : }
341 :
342 : /* Compute m-th derivative of inverse Mellin at t by continued fraction of
343 : * asymptotic expansion; t2d = t^(2/d). If t is NULL, "lfun" mode: don't
344 : * bother about complex branches + use absolute (rather than relative)
345 : * accuracy */
346 : static GEN
347 5620557 : Kderivlarge(GEN K, GEN t, GEN t2d, long bitprec0)
348 : {
349 : GEN tdA, P, S, pi, z;
350 5620557 : const long d = GMi_get_degree(K);
351 5620467 : GEN M = GMi_get_M(K), cd = GMi_get_cd(K), A2 = GMi_get_A2(K);
352 5620410 : long prec, nlim, status = GMi_get_status(K), m = GMi_get_m(K), bitprec = bitprec0;
353 :
354 5620365 : Kderivlarge_optim(K, !t, t2d, gtodouble(cd), &bitprec, &nlim);
355 5620105 : if (bitprec <= 0) return gen_0;
356 5534971 : prec = nbits2prec(bitprec);
357 5534946 : t2d = gtofp(t2d, prec);
358 5535452 : if (t)
359 4702 : tdA = gpow(t, gdivgu(A2,d), prec);
360 : else
361 5530750 : tdA = gpow(t2d, gdivgu(A2,2), prec);
362 5534987 : pi = mppi(prec); z = gmul(pi, t2d);
363 5535146 : P = gmul(gmul(cd, tdA), gexp(gmulsg(-d, z), prec));
364 5535297 : if (m) P = gmul(P, gpowgs(mulsr(-2, pi), m));
365 5535289 : if (status == 2) /* finite continued fraction */
366 1191932 : S = (lg(M) == 2)? gel(M,1): poleval(M, ginv(z));
367 : else
368 : {
369 4343357 : S = contfraceval_inv(M, z, nlim/2);
370 4343480 : if (DEBUGLEVEL>3)
371 : {
372 0 : GEN S0 = contfraceval_inv(M, z, minss(nlim/2 + 1, minss(lg(gel(M, 1)) - 1, lg(gel(M, 2)))));
373 0 : long e = gexpo(gmul(P, gsub(S,S0)));
374 0 : if (-e < bitprec0)
375 0 : err_printf("Kderivlarge: e = %ld, bit = %ld\n",e,bitprec0);
376 : }
377 4343480 : if (status == 1) S = gmul(S, gsubsg(1, ginv(gmul(z, pi))));
378 : }
379 5535411 : return gmul(P, S);
380 : }
381 :
382 : /* Dokchitser's coefficients used for asymptotic expansion of inverse Mellin
383 : * 2 <= p <= min(n+1, d), c = 2n-p+1; sh = (sh(x)/x)^(d-p) */
384 : static GEN
385 919182 : vp(long p, long c, GEN SMd, GEN sh)
386 : {
387 919182 : GEN s, ve = cgetg(p+2, t_VEC);
388 : long m, j, k;
389 :
390 919182 : gel(ve,1) = gen_1; gel(ve,2) = utoipos(c);
391 2690684 : for (j = 2; j <= p; j++) gel(ve,j+1) = gdivgs(gmulgu(gel(ve,j), c), j);
392 919182 : s = gel(SMd, 1);
393 3609866 : for (m = 1; m <= p; m++)
394 : {
395 2690684 : GEN t, c = gel(SMd, m+1);
396 2690684 : if (gequal0(c)) continue;
397 1409728 : t = gel(ve, m+1);
398 2825847 : for (k = 1; k <= m/2; k++)
399 1416119 : t = gadd(t, gmul(gel(ve, m-2*k+1), RgX_coeff(sh, k)));
400 1409728 : s = gadd(s, gmul(c, t));
401 : }
402 919182 : return s;
403 : }
404 :
405 : static GEN
406 5733 : get_SM(GEN Vga)
407 : {
408 5733 : long k, m, d = lg(Vga)-1;
409 5733 : GEN pol, nS1, SM, C, t = vecsum(Vga);
410 :
411 5733 : pol = roots_to_pol(gmulgs(Vga, -d), 0); /* deg(pol) = d */
412 5733 : SM = cgetg(d+2, t_VEC); gel(SM,1) = gen_1;
413 5733 : if (gequal0(t))
414 : { /* shortcut */
415 9688 : for (m = 1; m <= d; m++) gel(SM,m+1) = gel(pol,d+2-m);
416 2436 : return SM;
417 : }
418 3297 : nS1 = gpowers(gneg(t), d); C = matpascal(d);
419 14462 : for (m = 1; m <= d; m++)
420 : {
421 11165 : pari_sp av = avma;
422 11165 : GEN s = gmul(gel(nS1, m+1), gcoeff(C, d+1, m+1));
423 37660 : for (k = 1; k <= m; k++)
424 : {
425 26495 : GEN e = gmul(gel(nS1, m-k+1), gcoeff(C, d-k+1, m-k+1));
426 26495 : s = gadd(s, gmul(e, RgX_coeff(pol, d-k)));
427 : }
428 11165 : gel(SM, m+1) = gerepileupto(av, s);
429 : }
430 3297 : return SM;
431 : }
432 :
433 : static GEN
434 5733 : get_SMd(GEN Vga)
435 : {
436 5733 : GEN M, SM = get_SM(Vga);
437 5733 : long p, m, d = lg(Vga)-1;
438 :
439 5733 : M = cgetg(d, t_VEC);
440 18417 : for (p = 2; p <= d; p++)
441 : {
442 12684 : GEN a = gen_1, c;
443 12684 : long D = d - p;
444 12684 : gel(M, p-1) = c = cgetg(p+2, t_COL);
445 12684 : gel(c, 1) = gel(SM, p+1);
446 49504 : for (m = 1; m <= p; m++)
447 : {
448 36820 : a = muliu(a, D + m);
449 36820 : gel(c, m+1) = gmul(gel(SM, p-m+1), a);
450 : }
451 : }
452 5733 : return M;
453 : }
454 :
455 : /* Asymptotic expansion of inverse Mellin, to length nlimmax. Set status = 0
456 : * (regular), 1 (one Hankel determinant vanishes => contfracinit will fail)
457 : * or 2 (same as 1, but asymptotic expansion is finite!)
458 : *
459 : * If status = 2, the asymptotic expansion is finite so return only
460 : * the necessary number of terms nlim <= nlimmax + d. */
461 : static GEN
462 24157 : Klargeinit(GEN Vga, long nlimmax, long *status, long prec)
463 : {
464 24157 : long d = lg(Vga) - 1, p, n, cnt;
465 : GEN M, SMd, se, vsinh, vd;
466 :
467 24157 : if (Vgaeasytheta(Vga)) { *status = 2; return mkvec(gen_1); }
468 : /* d >= 2 */
469 5733 : *status = 0;
470 5733 : if (prec) prec += nbits2extraprec((prec >> 1) + BITS_IN_LONG);
471 5733 : SMd = get_SMd(Vga);
472 5733 : se = gsinh(RgX_to_ser(pol_x(0), d+2), 0); setvalser(se,0);
473 5733 : se = gdeflate(se, 0, 2); /* se(x^2) = sinh(x)/x */
474 5733 : vsinh = gpowers(se, d);
475 5733 : vd = gpowers(utoipos(2*d), d);
476 5733 : M = cgetg(nlimmax + d + 1, t_VEC); gel(M,1) = gen_1;
477 413362 : for (n = 2, cnt = 0; n <= nlimmax || cnt; n++)
478 : {
479 413362 : pari_sp av = avma;
480 413362 : long ld = minss(d, n);
481 413362 : GEN s = gen_0;
482 1332544 : for (p = 2; p <= ld; p++)
483 : {
484 919182 : GEN z = vp(p, 2*n-1-p, gel(SMd, p-1), gel(vsinh, d-p+1));
485 919182 : s = gadd(s, gmul(gdiv(z, gel(vd, p+1)), gel(M, n+1-p)));
486 : }
487 413362 : if (prec && !isinexact(s)) s = gtofp(s, prec);
488 413362 : gel(M,n) = s = gerepileupto(av, gdivgs(s, 1-n));
489 413362 : if (gequal0(s))
490 : {
491 56 : cnt++; *status = 1;
492 56 : if (cnt >= d-1) { *status = 2; n -= d-2; break; }
493 : }
494 : else
495 : {
496 413306 : if (n >= nlimmax) { n++; break; }
497 407601 : cnt = 0;
498 : }
499 : }
500 5733 : setlg(M, n); return M;
501 : }
502 :
503 : /* remove trailing zeros from vector. */
504 : static void
505 252 : stripzeros(GEN M)
506 : {
507 : long i;
508 441 : for(i = lg(M)-1; i >= 1; --i)
509 441 : if (!gequal0(gel(M, i))) break;
510 252 : setlg(M, i+1);
511 252 : }
512 :
513 : /* Asymptotic expansion of the m-th derivative of inverse Mellin, to length
514 : * nlimmax. If status = 2, the asymptotic expansion is finite so return only
515 : * the necessary number of terms nlim <= nlimmax + d. */
516 : static GEN
517 21371 : gammamellininvasymp_i(GEN Vga, long nlimmax, long m, long *status, long prec)
518 : {
519 : GEN M, A, Aadd;
520 : long d, i, nlim, n;
521 :
522 21371 : M = Klargeinit(Vga, nlimmax, status, prec);
523 21371 : if (!m) return M;
524 252 : d = lg(Vga)-1;
525 : /* half the exponent of t in asymptotic expansion. */
526 252 : A = gdivgu(gaddsg(1-d, sumVga(Vga)), 2*d);
527 252 : if (*status == 2) M = shallowconcat(M, zerovec(m));
528 252 : nlim = lg(M)-1;
529 252 : Aadd = sstoQ(2-d, 2*d); /* (1/d) - (1/2) */
530 532 : for (i = 1; i <= m; i++, A = gadd(A,Aadd))
531 8022 : for (n = nlim-1; n >= 1; --n)
532 15484 : gel(M, n+1) = gsub(gel(M, n+1),
533 15484 : gmul(gel(M, n), gsub(A, uutoQ(n-1, d))));
534 252 : stripzeros(M); return M;
535 : }
536 :
537 : INLINE int
538 21357 : RgV_is_CV(GEN x)
539 : {
540 : long i;
541 81655 : for (i = lg(x)-1; i > 0; i--)
542 : {
543 81067 : long t = typ(gel(x,i));
544 81067 : if (!is_real_t(t) && t!= t_COMPLEX) return 0;
545 : }
546 588 : return 1;
547 : }
548 :
549 : static GEN
550 21385 : get_Vga(GEN x, GEN *ldata)
551 : {
552 21385 : if (typ(x)==t_VEC && RgV_is_CV(x)) { *ldata = NULL; return x; }
553 20797 : *ldata = lfunmisc_to_ldata_shallow_i(x);
554 20797 : if (*ldata) x = ldata_get_gammavec(*ldata);
555 20797 : return x;
556 : }
557 : GEN
558 28 : gammamellininvasymp(GEN Vga, long nlim, long m)
559 : {
560 28 : pari_sp av = avma;
561 : long status;
562 : GEN ldata;
563 28 : Vga = get_Vga(Vga, &ldata);
564 28 : if (!is_vec_t(typ(Vga)) || lg(Vga) == 1)
565 7 : pari_err_TYPE("gammamellininvasymp",Vga);
566 21 : return gerepilecopy(av, gammamellininvasymp_i(Vga, nlim, m, &status, 0));
567 : }
568 :
569 : /* Does the continued fraction of the asymptotic expansion M at oo of inverse
570 : * Mellin transform attached to Vga have zero Hankel determinants ? */
571 : static long
572 2905 : ishankelspec(GEN Vga)
573 : {
574 2905 : long status, i, d = lg(Vga)-1;
575 : GEN M;
576 :
577 2905 : if (d == 5 || d == 7)
578 21 : { /* known bad cases: a x 5 and a x 7 */
579 133 : GEN v1 = gel(Vga, 1);
580 588 : for (i = 2; i <= d; ++i)
581 476 : if (!gequal(gel(Vga,i), v1)) break;
582 133 : if (i > d) return 1;
583 : }
584 2772 : else if (d==10 || d==14)
585 : { /* [ a x 5 , (a+1) x 5 ] */
586 7 : long d2 = d>>1;
587 7 : long s0 = 1, s1 = 0, sm1 = 0;
588 7 : GEN v1 = gel(Vga, 1);
589 70 : for (i = 2; i <= d; i++)
590 : {
591 63 : GEN s = gsub(gel(Vga,i),v1);
592 63 : if (gequal0(s)) s0++;
593 35 : else if(gequal1(s)) s1++;
594 0 : else if(gequalm1(s)) sm1++;
595 : }
596 7 : if (s0==d2 && (s1==d2 || sm1==d2)) return 1;
597 : }
598 : /* Heuristic: if 6 first terms in contfracinit don't fail, assume OK */
599 2786 : M = Klargeinit(Vga, 7, &status, 0);
600 2786 : return !contfracinit_i(M, 6);
601 : }
602 :
603 : /* Initialize data for computing m-th derivative of inverse Mellin */
604 : GEN
605 21364 : gammamellininvinit(GEN Vga, long m, long bitprec)
606 : {
607 21364 : const double C2 = MELLININV_CUTOFF;
608 21364 : pari_sp ltop = avma;
609 : GEN A2, M, VS, VL, cd, ldata;
610 21364 : long nlimmax, status, d, prec = nbits2prec((4*bitprec)/3);
611 21364 : double E = M_LN2*bitprec, tmax = get_tmax(bitprec); /* = E/C2 */
612 :
613 21364 : if (m < 0)
614 7 : pari_err_DOMAIN("gammamellininvinit", "derivation order", "<", gen_0, stoi(m));
615 21357 : Vga = get_Vga(Vga, &ldata); d = lg(Vga)-1;
616 21357 : if (!is_vec_t(typ(Vga)) || !d) pari_err_TYPE("gammamellininvinit",Vga);
617 21350 : nlimmax = ceil(E * log2(1+M_PI*tmax) * C2 / get_D(d));
618 21350 : A2 = gaddsg(m*(2-d) + 1-d, sumVga(Vga));
619 21350 : cd = (d <= 2)? gen_2: gsqrt(gdivgu(int2n(d+1), d), nbits2prec(bitprec));
620 : /* if in Klarge, we have |t| > tmax = E/C2, thus nlim < E*C2/D. */
621 21350 : M = gammamellininvasymp_i(Vga, nlimmax, m, &status, prec);
622 21350 : if (status == 2)
623 : {
624 18438 : tmax = -1.; /* only use Klarge */
625 18438 : VS = gen_0;
626 : }
627 : else
628 : {
629 2912 : VS = Kderivsmallinit(ldata, Vga, m, bitprec);
630 2912 : if (status == 0 && ishankelspec(Vga)) status = 1;
631 2912 : if (status == 1)
632 : { /* a Hankel determinant vanishes => contfracinit is undefined.
633 : So compute K(t) / (1 - 1/(pi^2*t)) instead of K(t)*/
634 133 : GEN t = ginv(mppi(prec));
635 : long i;
636 20902 : for (i = 2; i < lg(M); ++i)
637 20769 : gel(M, i) = gadd(gel(M, i), gmul(gel(M, i-1), t));
638 : }
639 : else
640 2779 : M = RgC_gtofp(M, prec); /* convert from rationals to t_REAL: faster */
641 2912 : M = contfracinit(M, lg(M)-2);
642 : }
643 21350 : VL = mkvec3(mkvec2(M, stoi(status)), cd, A2);
644 21350 : return gerepilecopy(ltop, mkvec5(dbltor(tmax), Vga, stoi(m), VS, VL));
645 : }
646 :
647 : /* Compute m-th derivative of inverse Mellin at s2d = s^(d/2) using
648 : * initialization data. Use Taylor expansion at 0 for |s2d| < tmax, and
649 : * asymptotic expansion at oo otherwise. WARNING: assume that accuracy
650 : * has been increased according to tmax by the CALLING program. */
651 : static GEN
652 5696193 : gammamellininvrt_i(GEN K, GEN s, GEN s2d, long bit)
653 : {
654 5696193 : if (dblmodulus(s2d) < GMi_get_tmax(K, bit))
655 75377 : return Kderivsmall(K, s, s2d, bit);
656 : else
657 5620503 : return Kderivlarge(K, s, s2d, bit);
658 : }
659 : GEN
660 5690362 : gammamellininvrt(GEN K, GEN s2d, long bit)
661 5690362 : { return gammamellininvrt_i(K, NULL, s2d, bit); }
662 :
663 : /* Compute inverse Mellin at s. K from gammamellininv OR a Vga, in which
664 : * case the initialization data is computed. */
665 : GEN
666 5824 : gammamellininv(GEN K, GEN s, long m, long bitprec)
667 : {
668 5824 : pari_sp av = avma;
669 : GEN s2d;
670 : long d;
671 :
672 5824 : if (!is_vec_t(typ(K)) || lg(K) != 6 || !is_vec_t(typ(GMi_get_Vga(K))))
673 98 : K = gammamellininvinit(K, m, bitprec);
674 5824 : d = GMi_get_degree(K);
675 5824 : s2d = gpow(s, gdivgu(gen_2, d), nbits2prec(bitprec));
676 5824 : return gerepileupto(av, gammamellininvrt_i(K, s, s2d, bitprec));
677 : }
|
