Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is to exceed 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - basemath - mellininv.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.12.1 lcov report (development 25685-2ac6d95b9c) Lines: 357 365 97.8 %
Date: 2020-08-06 06:06:08 Functions: 37 37 100.0 %
Legend: Lines: hit not hit

          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. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : 
      17             : /*******************************************************************/
      18             : /*               Computation of inverse Mellin                     */
      19             : /*               transforms of gamma products.                     */
      20             : /*******************************************************************/
      21             : #ifndef M_E
      22             : #define M_E 2.7182818284590452354
      23             : #endif
      24             : 
      25             : /* Handle complex Vga whose sum is real */
      26             : static GEN
      27       14364 : sumVga(GEN Vga) { return real_i(vecsum(Vga)); }
      28             : 
      29             : /* rough approximation to W0(a > -1/e), < 1% relative error */
      30             : double
      31       87863 : dbllambertW0(double a)
      32             : {
      33       87863 :   if (a < -0.2583)
      34             :   {
      35           0 :     const double c2 = -1./3, c3 = 11./72, c4 = -43./540, c5 = 769./17280;
      36           0 :     double p = sqrt(2*(M_E*a+1));
      37           0 :     if (a < -0.3243) return -1+p*(1+p*(c2+p*c3));
      38           0 :     return -1+p*(1+p*(c2+p*(c3+p*(c4+p*c5))));
      39             :   }
      40             :   else
      41             :   {
      42       87863 :     double Wd = log(1.+a);
      43       87863 :     Wd *= (1.-log(Wd/a))/(1.+Wd);
      44       87863 :     if (a < 0.6482 && a > -0.1838) return Wd;
      45       87863 :     return Wd*(1.-log(Wd/a))/(1.+Wd);
      46             :   }
      47             : }
      48             : 
      49             : /* rough approximation to W_{-1}(0 > a > -1/e), < 1% relative error */
      50             : double
      51       54250 : dbllambertW_1(double a)
      52             : {
      53       54250 :   if (a < -0.2464)
      54             :   {
      55         280 :     const double c2 = -1./3, c3 = 11./72, c4 = -43./540, c5 = 769./17280;
      56         280 :     double p = -sqrt(2*(M_E*a+1));
      57         280 :     if (a < -0.3243) return -1+p*(1+p*(c2+p*c3));
      58         175 :     return -1+p*(1+p*(c2+p*(c3+p*(c4+p*c5))));
      59             :   }
      60             :   else
      61             :   {
      62             :     double Wd;
      63       53970 :     a = -a; Wd = -log(a);
      64       53970 :     Wd *= (1.-log(Wd/a))/(1.-Wd);
      65       53970 :     if (a < 0.0056) return -Wd;
      66          21 :     return -Wd*(1.-log(Wd/a))/(1.-Wd);
      67             :   }
      68             : }
      69             : 
      70             : /* ac != 0 */
      71             : static double
      72      492919 : lemma526_i(double ac, double c, double t, double B)
      73             : {
      74      492919 :   double D = -B/ac; /* sgn(t) = sgn(a) = - sgn(D) */
      75      492919 :   if (D <= 0)
      76             :   {
      77      428743 :     if (D > -100)
      78             :     {
      79       54579 :       D = -exp(D) / t;
      80       54579 :       if (D < - 1/M_E) return 0;
      81       54250 :       D = dbllambertW_1(D);
      82             :     }
      83             :     else
      84             :     { /* avoid underflow, use asymptotic expansion */
      85      374164 :       double U = D - log(t);
      86      374164 :       D = U - log(-U);
      87             :     }
      88      428414 :     return pow(maxdd(t, -t * D), c);
      89             :   }
      90             :   else
      91             :   {
      92       64176 :     if (D < 100)
      93       10157 :       D = dbllambertW0(-exp(D) / t);
      94             :     else
      95             :     { /* avoid overflow, use asymptotic expansion */
      96       54019 :       double U = D - log(-t);
      97       54019 :       D = U - log(U);
      98             :     }
      99       64176 :     return pow(-t * D, c);
     100             :   }
     101             : }
     102             : /* b > 0, c > 0; solve x^a exp(-b x^(1/c)) < e^(-B) for x >= 0 */
     103             : double
     104          14 : dbllemma526(double a, double b, double c, double B)
     105             : {
     106             :   double ac;
     107          14 :   if (!a) return B <= 0? 0: pow(B/b, c);
     108          14 :   ac = a*c; if (B < 0) B = 1e-9;
     109          14 :   return lemma526_i(ac, c, ac/b, B);
     110             : }
     111             : /* Same, special case b/c = 2Pi, the only one needed: for c = d/2 */
     112             : double
     113     1654931 : dblcoro526(double a, double c, double B)
     114             : {
     115     1654931 :   if (!a) return B <= 0? 0: pow(B/(2*M_PI*c), c);
     116      492905 :   if (B < 0) B = 1e-9;
     117      492905 :   return lemma526_i(a*c, c, a/(2*M_PI), B);
     118             : }
     119             : 
     120             : static const double MELLININV_CUTOFF = 121.; /* C*C */
     121             : 
     122             : /* x real */
     123             : static GEN
     124        6034 : RMOD2(GEN x) { return gsub(x, gmul2n(gdiventgs(x,2), 1)); }
     125             : /* x real or complex, return canonical representative for x mod 2Z */
     126             : static GEN
     127        6034 : MOD2(GEN x)
     128        6034 : { return typ(x) == t_COMPLEX? mkcomplex(RMOD2(gel(x,1)), gel(x,2)): RMOD2(x); }
     129             : static GEN
     130        1946 : RgV_MOD2(GEN x)
     131        7980 : { pari_APPLY_same(MOD2(gel(x,i))); }
     132             : 
     133             : /* classes of poles of the gamma factor mod 2Z, sorted by increasing
     134             :  * Re(s) mod 2 (in [0,2[).*/
     135             : static GEN
     136        1946 : gammapoles(GEN Vga, long *pdV, long bit)
     137             : {
     138        1946 :   long i, m, emax, l = lg(Vga);
     139        1946 :   GEN P, B = RgV_MOD2(Vga), V = cgetg(l, t_VEC);
     140        1946 :   P = gen_indexsort(B, (void*)lexcmp, cmp_nodata);
     141        5103 :   for (i = m = 1; i < l;)
     142             :   {
     143        3157 :     GEN u = gel(B, P[i]);
     144             :     long k;
     145        6034 :     for(k = i+1; k < l; k++)
     146             :     {
     147        4088 :       GEN v = gsub(u, gel(B, P[k]));
     148        4088 :       if (!gequal0(v) && (!isinexactreal(v) || gexpo(v) > -bit)) break;
     149             :     }
     150        3157 :     gel(V, m++) = vecslice(P,i,k-1);
     151        3157 :     i = k;
     152             :   }
     153        1946 :   setlg(V, m); emax = 0;
     154        5103 :   for (i = 1; i < m; i++)
     155             :   {
     156        3157 :     long j, e = 0, li = lg(gel(V,i))-1;
     157        3157 :     GEN b = gel(B, gel(V,i)[1]);
     158        8890 :     for (j = 1; j < m; j++)
     159        5733 :       if (j != i) e -= gexpo(gsub(gel(B, gel(V,j)[1]), b));
     160        3157 :     emax = maxss(emax, e * li);
     161             :   }
     162        5103 :   for (i = 1; i < m; i++) gel(V,i) = vecpermute(Vga, gel(V,i));
     163        1946 :   *pdV = emax; return V;
     164             : }
     165             : 
     166             : static GEN
     167      193739 : sercoeff(GEN x, long n, long prec)
     168             : {
     169      193739 :   long N = n - valp(x);
     170      193739 :   return (N < 0)? gen_0: gprec_wtrunc(gel(x, N+2), prec);
     171             : }
     172             : 
     173             : /* prod_i Gamma(s/2 + (m+LA[i])/2), set t *= prod_i (s/2 + (m+LA[i])/2) */
     174             : static GEN
     175        5733 : get_gamma(GEN *pt, GEN LA, GEN m, int round, long precdl, long prec)
     176             : {
     177        5733 :   long i, l = lg(LA);
     178        5733 :   GEN pr = NULL, t = *pt;
     179       16142 :   for (i = 1; i < l; i++)
     180             :   {
     181       10409 :     GEN u, g, a = gmul2n(gadd(m, gel(LA,i)), -1);
     182       10409 :     if (round) a = ground(a);
     183       10409 :     u = deg1pol_shallow(ghalf, a, 0);
     184       10409 :     g = ggamma(RgX_to_ser(u, precdl), prec);
     185       10409 :     pr = pr? gmul(pr, g): g;
     186       10409 :     t = t? gmul(t, u): u;
     187             :   }
     188        5733 :   *pt = t; return pr;
     189             : }
     190             : /* generalized power series expansion of inverse Mellin around x = 0;
     191             :  * m-th derivative */
     192             : static GEN
     193        1946 : Kderivsmallinit(GEN ldata, GEN Vga, long m, long bit)
     194             : {
     195        1946 :   const double C2 = MELLININV_CUTOFF;
     196        1946 :   long prec2, N, j, l, dLA, limn, d = lg(Vga)-1;
     197             :   GEN piA, LA, L, M, mat;
     198             : 
     199        1946 :   LA = gammapoles(Vga, &dLA, bit); N = lg(LA)-1;
     200        1946 :   prec2 = nbits2prec(dLA + bit * (1 + M_PI*d/C2));
     201             : #if BITS_IN_LONG == 32
     202         278 :   prec2 += prec2 & 1L;
     203             : #endif
     204        1946 :   if (ldata) Vga = ldata_get_gammavec(ldata_newprec(ldata, prec2));
     205        1946 :   L = cgetg(N+1, t_VECSMALL);
     206        1946 :   M = cgetg(N+1, t_VEC);
     207        1946 :   mat = cgetg(N+1, t_VEC);
     208        1946 :   limn = ceil(2*M_LN2*bit / (d * dbllambertW0(C2/(M_PI*M_E))));
     209        1946 :   l = limn + 2;
     210        5103 :   for (j = 1; j <= N; j++)
     211             :   {
     212        3157 :     GEN S = gel(LA,j);
     213        3157 :     GEN C, c, mj, G = NULL, t = NULL, tj = NULL;
     214        3157 :     long i, k, n, jj, lj = L[j] = lg(S)-1, precdl = lj+3;
     215             : 
     216        3157 :     gel(M,j) = mj = gsubsg(2, gel(S, vecindexmin(real_i(S))));
     217        8890 :     for (jj = 1; jj <= N; jj++)
     218             :     {
     219             :       GEN g;
     220        5733 :       if (jj == j) /* poles come from this class only */
     221        3157 :         g = get_gamma(&tj, gel(LA,jj), mj, 1, precdl, prec2);
     222             :       else
     223        2576 :         g = get_gamma(&t, gel(LA,jj), mj, 0, precdl, prec2);
     224        5733 :       G = G? gmul(G, g): g;
     225             :     }
     226        3157 :     c = cgetg(limn+2,t_COL); gel(c,1) = G;
     227      108738 :     for (n=1; n <= limn; n++)
     228             :     {
     229      105581 :       GEN A = utoineg(2*n), T = RgX_translate(tj, A);
     230             :       /* T = exact polynomial, may vanish at 0 (=> pole in c[n+1]) */
     231      105581 :       if (t) T = RgX_mul(T, RgX_translate(t, A)); /* no pole here */
     232      105581 :       gel(c,n+1) = gdiv(gel(c,n), T);
     233             :     }
     234        3157 :     gel(mat, j) = C = cgetg(lj+1, t_COL);
     235        9191 :     for (k = 1; k <= lj; k++)
     236             :     {
     237        6034 :       GEN L = cgetg(l, t_POL);
     238      199773 :       for (n = 2; n < l; n++) gel(L,n) = sercoeff(gel(c,n), -k, prec2);
     239        6034 :       L[1] = evalsigne(1)|evalvarn(0); gel(C,k) = L;
     240             :     }
     241             :     /* C[k] = \sum_n c_{j,k} t^n =: C_k(t) in Dokchitser's Algo 3.3
     242             :      * m-th derivative of t^(-M+2) sum_k (-ln t)^k/k! C_k(t^2) */
     243        3157 :     if (m)
     244             :     {
     245         119 :       mj = gsubgs(mj, 2);
     246         238 :       for (i = 1; i <= m; i++, mj = gaddgs(mj,1))
     247         322 :         for (k = 1; k <= lj; k++)
     248             :         {
     249         203 :           GEN c = gel(C,k), d = RgX_shift_shallow(gmul2n(RgX_deriv(c),1), 1);
     250         203 :           c = RgX_Rg_mul(c, mj);
     251         203 :           if (k < lj) c = RgX_add(c, gel(C,k+1));
     252         203 :           gel(C,k) = RgX_sub(d, c);
     253             :         }
     254         119 :       gel(M,j) = gaddgs(mj,2);
     255             :     }
     256        9191 :     for (k = 1; k <= lj; k++)
     257             :     {
     258        6034 :       GEN c = gel(C,k);
     259        6034 :       if (k > 2) c = RgX_Rg_div(c, mpfact(k-1));
     260        6034 :       gel(C,k) = RgX_to_RgC(c, lgpol(c));
     261             :     }
     262             :   }
     263             :   /* Algo 3.3: * \phi^(m)(t) = sum_j t^m_j sum_k (-ln t)^k mat[j,k](t^2) */
     264        1946 :   piA = gsubsg(m*d, sumVga(Vga));
     265        1946 :   if (!gequal0(piA)) piA = powPis(gmul2n(piA,-1), prec2);
     266        1946 :   return mkvec5(L, RgV_neg(M), mat, mkvecsmall(prec2), piA);
     267             : }
     268             : 
     269             : /* Evaluate a vector considered as a polynomial using Horner. Unstable!
     270             :  * If ui != NULL, ui = 1/u, evaluate P(1/u)*u^(deg P): useful for |u|>1 */
     271             : static GEN
     272      178469 : evalvec(GEN vec, long lim, GEN u, GEN ui)
     273             : {
     274      178469 :   pari_sp ltop = avma;
     275      178469 :   GEN S = gen_0;
     276             :   long n;
     277      178469 :   lim = minss(lim, lg(vec)-1);
     278      178469 :   if (!ui)
     279     1281019 :     for (n = lim; n >= 1; --n) S = gmul(u, gadd(gel(vec,n), S));
     280             :   else
     281             :   {
     282     3984957 :     for (n = 1; n <= lim; ++n) S = gmul(ui, gadd(gel(vec,n), S));
     283      126623 :     S = gmul(gpowgs(u, n), S);
     284             :   }
     285      178466 :   return gerepileupto(ltop, S);
     286             : }
     287             : 
     288             : /* gammamellininvinit accessors */
     289             : static double
     290     2963509 : get_tmax(long bitprec)
     291     2963509 : { return (M_LN2 / MELLININV_CUTOFF) * bitprec ; }
     292             : static GEN
     293        5726 : GMi_get_Vga(GEN K) { return gel(K,2); }
     294             : static long
     295     6644142 : GMi_get_degree(GEN K) { return lg(gel(K,2))-1; }
     296             : static long
     297     3293795 : GMi_get_m(GEN K) { return itos( gel(K,3) ); }
     298             : static GEN /* [lj,mj,mat,prec2], Kderivsmall only */
     299     3397034 : GMi_get_VS(GEN K) { return gel(K,4); }
     300             : static GEN /* [Ms,cd,A2], Kderivlarge only */
     301     6586641 : GMi_get_VL(GEN K) { return gel(K,5); }
     302             : static double
     303     3345172 : GMi_get_tmax(GEN K, long bitprec)
     304     3345172 : { return (typ(GMi_get_VS(K)) == t_INT)? -1.0 : get_tmax(bitprec); }
     305             : 
     306             : /* Compute m-th derivative of inverse Mellin at x by generalized power series
     307             :  * around x = 0; x2d = x^(2/d), x is possibly NULL (don't bother about
     308             :  * complex branches). Assume |x|^(2/d) <= tmax = M_LN2*bitprec/MELLININV_CUTOFF*/
     309             : static GEN
     310       51862 : Kderivsmall(GEN K, GEN x, GEN x2d, long bitprec)
     311             : {
     312       51862 :   GEN VS = GMi_get_VS(K), L = gel(VS,1), M = gel(VS,2), mat = gel(VS,3);
     313       51862 :   GEN d2, Lx, x2, x2i, S, pi, piA = gel(VS,5);
     314       51862 :   long prec = gel(VS,4)[1], d = GMi_get_degree(K);
     315       51862 :   long j, k, limn, N = lg(L)-1;
     316       51862 :   double xd, Wd, Ed = M_LN2*bitprec / d;
     317             : 
     318       51862 :   xd = maxdd(M_PI*dblmodulus(x2d), 1E-13); /* pi |x|^2/d unless x tiny */
     319       51862 :   if (xd > Ed) pari_err_BUG("Kderivsmall (x2d too large)");
     320             :   /* Lemma 5.2.6 (2), a = 1 + log(Pi x^(2/d)) = log(e / xd),
     321             :    * B = log(2)*bitprec / d = Ed */
     322       51862 :   Wd = dbllambertW0( Ed / (M_E*xd) ); /* solution of w exp(w) = B exp(-a)*/
     323       51862 :   limn = (long) ceil(2*Ed/Wd);
     324       51862 :   pi = mppi(prec);
     325       51862 :   d2 = gdivsg(d,gen_2);
     326       51860 :   if (x)
     327        1122 :     x = gmul(gtofp(x,prec), gpow(pi,d2,prec));
     328             :   else
     329       50738 :     x = gpow(gmul(gtofp(x2d,prec),pi), d2, prec);
     330             :   /* at this stage, x has been replaced by pi^(d/2) x */
     331       51862 :   x2 = gsqr(x);
     332       51863 :   Lx = gpowers(gneg(glog(x,prec)), vecsmall_max(L));
     333       51862 :   x2i = (gcmp(gnorml2(x2), gen_1) <= 0)? NULL: ginv(x2);
     334       51860 :   S = gen_0;
     335      134803 :   for (j = 1; j <= N; ++j)
     336             :   {
     337       82944 :     long lj = L[j];
     338       82944 :     GEN s = gen_0;
     339      261413 :     for (k = 1; k <= lj; k++)
     340      178469 :       s = gadd(s, gmul(gel(Lx,k), evalvec(gmael(mat,j,k), limn, x2, x2i)));
     341       82944 :     S = gadd(S, gmul(gpow(x, gel(M,j), prec), s));
     342             :   }
     343       51859 :   if (!gequal0(piA)) S = gmul(S, piA);
     344       51859 :   return S;
     345             : }
     346             : 
     347             : /* In Klarge, we conpute K(t) as (asymptotic) * F(z), where F ~ 1 is given by
     348             :  * a continued fraction and z = Pi t^(2/d). If we take 2n terms in F (n terms
     349             :  * in Euler form), F_n(z) - F(z) is experimentally in exp(- C sqrt(n*z))
     350             :  * where C ~ 8 for d > 2 [HEURISTIC] and C = 4 (theorem) for d = 1 or d = 2
     351             :  * and vga = [0,1]. For e^(-E) absolute error, we want
     352             :  *   exp(-C sqrt(nz)) < e^-(E+a), where a ~ ln(asymptotic)
     353             :  * i.e.  2n > (E+a)^2 / t^(2/d) * 2/(C^2 Pi); C^2*Pi/2 ~ 100.5 ~ 101
     354             :  *
     355             :  * In fact, this model becomes wrong for z large: we use instead
     356             :  *
     357             :  *   exp(- sqrt(D * nz/log(z+1))) < e^-(E+a),
     358             :  * i.e.  2n > (E+a)^2 * log(1 + Pi t^(2/d))/ t^(2/d) * 2/(D Pi); */
     359             : static double
     360     3305885 : get_D(long d) { return d <= 2 ? 157. : 180.; }
     361             : /* if (abs), absolute error rather than relative */
     362             : static void
     363     3293769 : Kderivlarge_optim(GEN K, long abs, GEN t2d,GEN gcd, long *pbitprec, long *pnlim)
     364             : {
     365     3293769 :   GEN VL = GMi_get_VL(K), A2 = gel(VL,3);
     366     3293682 :   long bitprec = *pbitprec, d = GMi_get_degree(K);
     367     3293731 :   const double D = get_D(d), td = dblmodulus(t2d), cd = gtodouble(gcd);
     368     3293229 :   double a, rtd, E = M_LN2*bitprec;
     369             : 
     370     3293229 :   rtd = (typ(t2d) == t_COMPLEX)? gtodouble(gel(t2d,1)): td;
     371             : 
     372             :   /* A2/2 = A, log(td) = (2/d)*log t */
     373     3293229 :   a = d*gtodouble(A2)*log2(td)/2 - (M_PI/M_LN2)*d*rtd + log2(cd); /*log2 K(t)~a*/
     374             : 
     375             :   /* if bitprec <= 0, caller should return K(t) ~ 0 */
     376     3292904 :   bitprec += 64;
     377     3292904 :   if (abs)
     378             :   {
     379     3288337 :     bitprec += ceil(a);
     380     3288337 :     if (a <= -65) E = M_LN2*bitprec; /* guarantees E <= initial E */
     381             :   }
     382     3292904 :   *pbitprec = bitprec;
     383     3292904 :   *pnlim = ceil(E*E * log2(1+M_PI*td) / (D*td));
     384     3292904 : }
     385             : 
     386             : /* Compute m-th derivative of inverse Mellin at t by continued fraction of
     387             :  * asymptotic expansion; t2d = t^(2/d). If t is NULL, "lfun" mode: don't
     388             :  * bother about complex branches + use absolute (rather than relative)
     389             :  * accuracy */
     390             : static GEN
     391     3293796 : Kderivlarge(GEN K, GEN t, GEN t2d, long bitprec0)
     392             : {
     393             :   GEN tdA, P, S, pi, z;
     394     3293796 :   const long d = GMi_get_degree(K);
     395     3293621 :   GEN M, VL = GMi_get_VL(K), Ms = gel(VL,1), cd = gel(VL,2), A2 = gel(VL,3);
     396     3293594 :   long status, prec, nlim, m = GMi_get_m(K), bitprec = bitprec0;
     397             : 
     398     3293744 :   Kderivlarge_optim(K, !t, t2d, cd, &bitprec, &nlim);
     399     3293378 :   if (bitprec <= 0) return gen_0;
     400     3226185 :   prec = nbits2prec(bitprec);
     401     3226102 :   t2d = gtofp(t2d, prec);
     402     3226229 :   if (t)
     403        4695 :     tdA = gpow(t, gdivgs(A2,d), prec);
     404             :   else
     405     3221534 :     tdA = gpow(t2d, gdivgs(A2,2), prec);
     406     3225696 :   tdA = gmul(cd, tdA);
     407             : 
     408     3225311 :   pi = mppi(prec);
     409     3225592 :   z = gmul(pi, t2d);
     410     3226327 :   P = gmul(tdA, gexp(gmulsg(-d, z), prec));
     411     3226238 :   if (m) P = gmul(P, gpowgs(mulsr(-2, pi), m));
     412     3226238 :   M = gel(Ms,1);
     413     3226238 :   status = itos(gel(Ms,2));
     414     3226055 :   if (status == 2)
     415             :   {
     416      391542 :     if (lg(M) == 2) /* shortcut: continued fraction is constant */
     417      391213 :       S = gel(M,1);
     418             :     else
     419         329 :       S = poleval(RgV_to_RgX(M, 0), ginv(z));
     420             :   }
     421             :   else
     422             :   {
     423     2834513 :     S = contfraceval_inv(M, z, nlim/2);
     424     2834895 :     if (DEBUGLEVEL>3)
     425             :     {
     426           0 :       GEN S0 = contfraceval_inv(M, z, nlim/2 + 1);
     427           0 :       long e = gexpo(gmul(P, gabs(gsub(S,S0),0)));
     428           0 :       if (-e < bitprec0)
     429           0 :         err_printf("Kderivlarge: e = %ld, bit = %ld\n",e,bitprec0);
     430             :     }
     431     2834895 :     if (status == 1) S = gmul(S, gsubsg(1, ginv(gmul(z, pi))));
     432             :   }
     433     3226427 :   return gmul(P, S);
     434             : }
     435             : 
     436             : /* Dokchitser's coefficients used for asymptotic expansion of inverse Mellin
     437             :  * 2 <= p <= min(n+1, d), c = 2n-p+1; sh = (sh(x)/x)^(d-p) */
     438             : static GEN
     439      515312 : vp(long p, long c, GEN SMd, GEN sh)
     440             : {
     441      515312 :   GEN s, ve = cgetg(p+2, t_VEC);
     442             :   long m, j, k;
     443             : 
     444      515312 :   gel(ve,1) = gen_1; gel(ve,2) = utoipos(c);
     445     1425683 :   for (j = 2; j <= p; j++) gel(ve,j+1) = gdivgs(gmulgs(gel(ve,j), c), j);
     446      515312 :   s = gel(SMd, 1);
     447     1940995 :   for (m = 1; m <= p; m++)
     448             :   {
     449     1425683 :     GEN t, c = gel(SMd, m+1);
     450     1425683 :     if (gequal0(c)) continue;
     451      760153 :     t = gel(ve, m+1);
     452     1485507 :     for (k = 1; k <= m/2; k++)
     453      725354 :       t = gadd(t, gmul(gel(ve, m-2*k+1), RgX_coeff(sh, k)));
     454      760153 :     s = gadd(s, gmul(c, t));
     455             :   }
     456      515312 :   return s;
     457             : }
     458             : 
     459             : static GEN
     460        1981 : get_SM(GEN Vga)
     461             : {
     462        1981 :   long k, m, d = lg(Vga)-1;
     463        1981 :   GEN pol, nS1, SM, C, t = vecsum(Vga);
     464             : 
     465        1981 :   pol = roots_to_pol(gmulgs(Vga, -d), 0); /* deg(pol) = d */
     466        1981 :   SM = cgetg(d+2, t_VEC); gel(SM,1) = gen_1;
     467        1981 :   if (gequal0(t))
     468             :   { /* shortcut */
     469        3465 :     for (m = 1; m <= d; m++) gel(SM,m+1) = gel(pol,d+2-m);
     470         903 :     return SM;
     471             :   }
     472        1078 :   nS1 = gpowers(gneg(t), d); C = matpascal(d);
     473        4634 :   for (m = 1; m <= d; m++)
     474             :   {
     475        3556 :     pari_sp av = avma;
     476        3556 :     GEN s = gmul(gel(nS1, m+1), gcoeff(C, d+1, m+1));
     477       11781 :     for (k = 1; k <= m; k++)
     478             :     {
     479        8225 :       GEN e = gmul(gel(nS1, m-k+1), gcoeff(C, d-k+1, m-k+1));
     480        8225 :       s = gadd(s, gmul(e, RgX_coeff(pol, d-k)));
     481             :     }
     482        3556 :     gel(SM, m+1) = gerepileupto(av, s);
     483             :   }
     484        1078 :   return SM;
     485             : }
     486             : 
     487             : static GEN
     488        1981 : get_SMd(GEN Vga)
     489             : {
     490        1981 :   GEN M, SM = get_SM(Vga);
     491        1981 :   long p, m, d = lg(Vga)-1;
     492             : 
     493        1981 :   M = cgetg(d, t_VEC);
     494        6118 :   for (p = 2; p <= d; p++)
     495             :   {
     496        4137 :     GEN a = gen_1, c;
     497        4137 :     long D = d - p;
     498        4137 :     gel(M, p-1) = c = cgetg(p+2, t_COL);
     499        4137 :     gel(c, 1) = gel(SM, p+1);
     500       15785 :     for (m = 1; m <= p; m++)
     501             :     {
     502       11648 :       a = muliu(a, D + m);
     503       11648 :       gel(c, m+1) = gmul(gel(SM, p-m+1), a);
     504             :     }
     505             :   }
     506        1981 :   return M;
     507             : }
     508             : 
     509             : /* Asymptotic expansion of inverse Mellin, to length nlimmax. Set status = 0
     510             :  * (regular), 1 (one Hankel determinant vanishes => contfracinit will fail)
     511             :  * or 2 (same as 1, but asymptotic expansion is finite!)
     512             :  *
     513             :  * If status = 2, the asymptotic expansion is finite so return only
     514             :  * the necessary number of terms nlim <= nlimmax + d. */
     515             : static GEN
     516       12187 : Klargeinit(GEN Vga, long nlimmax, long *status)
     517             : {
     518       12187 :   long d = lg(Vga) - 1, p, n, cnt;
     519             :   GEN M, SMd, se, vsinh, vd;
     520             : 
     521       12187 :   if (Vgaeasytheta(Vga)) { *status = 2; return mkvec(gen_1); }
     522             :   /* d >= 2 */
     523        1981 :   *status = 0;
     524        1981 :   SMd = get_SMd(Vga);
     525        1981 :   se = gsinh(RgX_to_ser(pol_x(0), d+2), 0); setvalp(se,0);
     526        1981 :   se = gdeflate(se, 0, 2); /* se(x^2) = sinh(x)/x */
     527        1981 :   vsinh = gpowers(se, d);
     528        1981 :   vd = gpowers(utoipos(2*d), d);
     529        1981 :   M = cgetg(nlimmax + d + 1, t_VEC); gel(M,1) = gen_1;
     530      251846 :   for (n = 2, cnt = 0; n <= nlimmax || cnt; n++)
     531             :   {
     532      251846 :     pari_sp av = avma;
     533      251846 :     long ld = minss(d, n);
     534      251846 :     GEN s = gen_0;
     535      767158 :     for (p = 2; p <= ld; p++)
     536             :     {
     537      515312 :       GEN z = vp(p, 2*n-1-p, gel(SMd, p-1), gel(vsinh, d-p+1));
     538      515312 :       s = gadd(s, gmul(gdiv(z, gel(vd, p+1)), gel(M, n+1-p)));
     539             :     }
     540      251846 :     gel(M,n) = s = gerepileupto(av, gdivgs(s, 1-n));
     541      251846 :     if (gequal0(s))
     542             :     {
     543          49 :       cnt++; *status = 1;
     544          49 :       if (cnt >= d-1) { *status = 2; n -= d-2; break; }
     545             :     }
     546             :     else
     547             :     {
     548      251797 :       if (n >= nlimmax) { n++; break; }
     549      249844 :       cnt = 0;
     550             :     }
     551             :   }
     552        1981 :   setlg(M, n); return M;
     553             : }
     554             : 
     555             : /* remove trailing zeros from vector. */
     556             : static void
     557         252 : stripzeros(GEN M)
     558             : {
     559             :   long i;
     560         441 :   for(i = lg(M)-1; i >= 1; --i)
     561         441 :     if (!gequal0(gel(M, i))) break;
     562         252 :   setlg(M, i+1);
     563         252 : }
     564             : 
     565             : /* Asymptotic expansion of the m-th derivative of inverse Mellin, to length
     566             :  * nlimmax. If status = 2, the asymptotic expansion is finite so return only
     567             :  * the necessary number of terms nlim <= nlimmax + d. */
     568             : static GEN
     569       12187 : gammamellininvasymp_i(GEN Vga, long nlimmax, long m, long *status)
     570             : {
     571             :   GEN M, A, Aadd;
     572             :   long d, i, nlim, n;
     573             : 
     574       12187 :   M = Klargeinit(Vga, nlimmax, status);
     575       12187 :   if (!m) return M;
     576         252 :   d = lg(Vga)-1;
     577             :   /* half the exponent of t in asymptotic expansion. */
     578         252 :   A = gdivgs(gaddsg(1-d, sumVga(Vga)), 2*d);
     579         252 :   if (*status == 2) M = shallowconcat(M, zerovec(m));
     580         252 :   nlim = lg(M)-1;
     581         252 :   Aadd = sstoQ(2-d, 2*d); /* (1/d) - (1/2) */
     582         532 :   for (i = 1; i <= m; i++, A = gadd(A,Aadd))
     583        8022 :     for (n = nlim-1; n >= 1; --n)
     584       15484 :       gel(M, n+1) = gsub(gel(M, n+1),
     585        7742 :                          gmul(gel(M, n), gsub(A, sstoQ(n-1, d))));
     586         252 :   stripzeros(M); return M;
     587             : }
     588             : static GEN
     589       12201 : get_Vga(GEN x, GEN *ldata)
     590             : {
     591       12201 :   *ldata = lfunmisc_to_ldata_shallow_i(x);
     592       12201 :   if (*ldata) x = ldata_get_gammavec(*ldata);
     593       12201 :   return x;
     594             : }
     595             : GEN
     596          28 : gammamellininvasymp(GEN Vga, long nlim, long m)
     597             : {
     598          28 :   pari_sp av = avma;
     599             :   long status;
     600             :   GEN ldata;
     601          28 :   Vga = get_Vga(Vga, &ldata);
     602          28 :   if (!is_vec_t(typ(Vga)) || lg(Vga) == 1)
     603           7 :     pari_err_TYPE("gammamellininvasymp",Vga);
     604          21 :   return gerepilecopy(av, gammamellininvasymp_i(Vga, nlim, m, &status));
     605             : }
     606             : 
     607             : /* Does the continued fraction of the asymptotic expansion M at oo of inverse
     608             :  * Mellin transform attached to Vga have zero Hankel determinants ? */
     609             : static long
     610        1939 : ishankelspec(GEN Vga, GEN M)
     611             : {
     612        1939 :   long status, i, d = lg(Vga)-1;
     613             : 
     614        1939 :   if (d == 5 || d == 7)
     615             :   { /* known bad cases: a x 5 and a x 7 */
     616          70 :     GEN v1 = gel(Vga, 1);
     617         294 :     for (i = 2; i <= d; ++i)
     618         238 :       if (!gequal(gel(Vga,i), v1)) break;
     619          70 :     if (i > d) return 1;
     620             :   }
     621        1883 :   status = 0;
     622             :   /* Heuristic: if 6 first terms in contfracinit don't fail, assume it's OK */
     623        1883 :   pari_CATCH(e_INV) { status = 1; }
     624        1883 :   pari_TRY { contfracinit(M, minss(lg(M)-2,6)); }
     625        1883 :   pari_ENDCATCH; return status;
     626             : }
     627             : 
     628             : /* Initialize data for computing m-th derivative of inverse Mellin */
     629             : GEN
     630       12173 : gammamellininvinit(GEN Vga, long m, long bitprec)
     631             : {
     632       12173 :   const double C2 = MELLININV_CUTOFF;
     633       12173 :   pari_sp ltop = avma;
     634             :   GEN A2, M, VS, VL, cd, ldata;
     635       12173 :   long nlimmax, status, d, prec = nbits2prec((4*bitprec)/3);
     636       12173 :   double E = M_LN2*bitprec, tmax = get_tmax(bitprec); /* = E/C2 */
     637             : 
     638       12173 :   Vga = get_Vga(Vga, &ldata); d = lg(Vga)-1;
     639       12173 :   if (!is_vec_t(typ(Vga)) || !d) pari_err_TYPE("gammamellininvinit",Vga);
     640       12166 :   nlimmax = ceil(E * log2(1+M_PI*tmax) * C2 / get_D(d));
     641       12166 :   A2 = gaddsg(m*(2-d) + 1-d, sumVga(Vga));
     642       12166 :   cd = (d <= 2)? gen_2: gsqrt(gdivgs(int2n(d+1), d), nbits2prec(bitprec));
     643             :   /* if in Klarge, we have |t| > tmax = E/C2, thus nlim < E*C2/D. */
     644       12166 :   M = gammamellininvasymp_i(RgV_gtofp(Vga, prec), nlimmax, m, &status);
     645       12166 :   if (status == 2)
     646             :   {
     647       10220 :     tmax = -1.; /* only use Klarge */
     648       10220 :     VS = gen_0;
     649             :   }
     650             :   else
     651             :   {
     652        1946 :     VS = Kderivsmallinit(ldata, Vga, m, bitprec);
     653        1946 :     if (status == 0 && ishankelspec(Vga, M)) status = 1;
     654        1946 :     if (status == 1)
     655             :     { /* a Hankel determinant vanishes => contfracinit is undefined.
     656             :          So compute K(t) / (1 - 1/(pi^2*t)) instead of K(t)*/
     657          63 :       GEN t = ginv(mppi(prec));
     658             :       long i;
     659        9772 :       for (i = 2; i < lg(M); ++i)
     660        9709 :         gel(M, i) = gadd(gel(M, i), gmul(gel(M, i-1), t));
     661             :     }
     662             :     else
     663        1883 :       M = RgC_gtofp(M, prec); /* convert from rationals to t_REAL: faster */
     664        1946 :     M = contfracinit(M, lg(M)-2);
     665             :   }
     666       12166 :   VL = mkvec3(mkvec2(M, stoi(status)), cd, A2);
     667       12166 :   return gerepilecopy(ltop, mkvec5(dbltor(tmax), Vga, stoi(m), VS, VL));
     668             : }
     669             : 
     670             : /* Compute m-th derivative of inverse Mellin at s2d = s^(d/2) using
     671             :  * initialization data. Use Taylor expansion at 0 for |s2d| < tmax, and
     672             :  * asymptotic expansion at oo otherwise. WARNING: assume that accuracy
     673             :  * has been increased according to tmax by the CALLING program. */
     674             : GEN
     675     3339800 : gammamellininvrt(GEN K, GEN s2d, long bitprec)
     676             : {
     677     3339800 :   if (dblmodulus(s2d) < GMi_get_tmax(K, bitprec))
     678       50739 :     return Kderivsmall(K, NULL, s2d, bitprec);
     679             :   else
     680     3288607 :     return Kderivlarge(K, NULL, s2d, bitprec);
     681             : }
     682             : 
     683             : /* Compute inverse Mellin at s. K from gammamellininv OR a Vga, in which
     684             :  * case the initialization data is computed. */
     685             : GEN
     686        5817 : gammamellininv(GEN K, GEN s, long m, long bitprec)
     687             : {
     688        5817 :   pari_sp av = avma;
     689             :   GEN z, s2d;
     690             :   long d;
     691             : 
     692        5817 :   if (!is_vec_t(typ(K)) || lg(K) != 6 || !is_vec_t(typ(GMi_get_Vga(K))))
     693          91 :     K = gammamellininvinit(K, m, bitprec);
     694        5817 :   d = GMi_get_degree(K);
     695        5817 :   s2d = gpow(s, gdivgs(gen_2, d), nbits2prec(bitprec));
     696        5817 :   if (dblmodulus(s2d) < GMi_get_tmax(K, bitprec))
     697        1122 :     z = Kderivsmall(K, s, s2d, bitprec);
     698             :   else
     699        4695 :     z = Kderivlarge(K, s, s2d, bitprec);
     700        5817 :   return gerepileupto(av, z);
     701             : }

Generated by: LCOV version 1.13