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.0 lcov report (development 23011-59c7027a2) Lines: 315 327 96.3 %
Date: 2018-09-22 05:37:52 Functions: 30 31 96.8 %
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       71809 : sumVga(GEN Vga) { return real_i(vecsum(Vga)); }
      28             : 
      29             : /* rough approximation to W0(a > -1/e), < 1% relative error */
      30             : double
      31       70538 : dbllambertW0(double a)
      32             : {
      33       70538 :   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       70538 :     double Wd = log(1.+a);
      43       70538 :     Wd *= (1.-log(Wd/a))/(1.+Wd);
      44       70538 :     if (a < 0.6482 && a > -0.1838) return Wd;
      45       70538 :     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       53102 : dbllambertW_1(double a)
      52             : {
      53       53102 :   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       52822 :     a = -a; Wd = -log(a);
      64       52822 :     Wd *= (1.-log(Wd/a))/(1.-Wd);
      65       52822 :     if (a < 0.0056) return -Wd;
      66          42 :     return -Wd*(1.-log(Wd/a))/(1.-Wd);
      67             :   }
      68             : }
      69             : 
      70             : /* ac != 0 */
      71             : static double
      72      501808 : lemma526_i(double ac, double c, double t, double B)
      73             : {
      74      501808 :   double D = -B/ac; /* sgn(t) = sgn(a) = - sgn(D) */
      75      501808 :   if (D <= 0)
      76             :   {
      77      448350 :     if (D > -100)
      78             :     {
      79       53431 :       D = -exp(D) / t;
      80       53431 :       if (D < - 1/M_E) return 0;
      81       53102 :       D = dbllambertW_1(D);
      82             :     }
      83             :     else
      84             :     { /* avoid underflow, use asymptotic expansion */
      85      394919 :       double U = D - log(t);
      86      394919 :       D = U - log(-U);
      87             :     }
      88      448021 :     return pow(maxdd(t, -t * D), c);
      89             :   }
      90             :   else
      91             :   {
      92       53458 :     if (D < 100)
      93        8935 :       D = dbllambertW0(-exp(D) / t);
      94             :     else
      95             :     { /* avoid overflow, use asymptotic expansion */
      96       44523 :       double U = D - log(-t);
      97       44523 :       D = U - log(U);
      98             :     }
      99       53458 :     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           0 : dbllemma526(double a, double b, double c, double B)
     105             : {
     106             :   double ac;
     107           0 :   if (!a) return B <= 0? 0: pow(B/b, c);
     108           0 :   ac = a*c; if (B < 0) B = 1e-9;
     109           0 :   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     1517900 : dblcoro526(double a, double c, double B)
     114             : {
     115     1517900 :   if (!a) return B <= 0? 0: pow(B/(2*M_PI*c), c);
     116      501808 :   if (B < 0) B = 1e-9;
     117      501808 :   return lemma526_i(a*c, c, a/(2*M_PI), B);
     118             : }
     119             : 
     120             : static const double MELLININV_CUTOFF = 121.; /* C*C */
     121             : 
     122             : static GEN
     123        5565 : MOD2(GEN x) { GEN q = gdivent(real_i(x),gen_2); return gsub(x,gmul2n(q,1)); }
     124             : static GEN
     125        1827 : RgV_MOD2(GEN v)
     126             : {
     127             :   long i, l;
     128        1827 :   GEN w = cgetg_copy(v,&l);
     129        1827 :   for (i=1; i<l; i++) gel(w,i) = MOD2(gel(v,i));
     130        1827 :   return w;
     131             : }
     132             : 
     133             : /* poles of the gamma factor */
     134             : static GEN
     135        1827 : gammapoles(GEN Vga)
     136             : {
     137        1827 :   long i, m, d = lg(Vga)-1;
     138        1827 :   GEN P, V, B = RgV_MOD2(Vga);
     139        1827 :   P = gen_indexsort(real_i(B), (void*)gcmp, cmp_nodata);
     140        1827 :   V = cgetg(d+1, t_VEC);
     141        6496 :   for (i = 1, m = 1; i <= d;)
     142             :   {
     143        2842 :     GEN u = gel(B, P[i]);
     144             :     long k;
     145        5565 :     for(k = i+1; k <= d; ++k)
     146        3738 :       if (!gequal(gel(B, P[k]), u)) break;
     147        2842 :     gel(V, m++) = vecpermute(Vga, vecslice(P,i,k-1));
     148        2842 :     i = k;
     149             :   }
     150        1827 :   setlg(V, m); return V;
     151             : }
     152             : 
     153             : static GEN
     154      186345 : sercoeff(GEN x, long n, long prec)
     155             : {
     156      186345 :   long N = n - valp(x);
     157      186345 :   return (N < 0)? gen_0: gtofp(gel(x, N+2), prec);
     158             : }
     159             : 
     160             : /* generalized power series expansion of inverse Mellin around x = 0;
     161             :  * m-th derivative */
     162             : static GEN
     163        1827 : Kderivsmallinit(GEN Vga, long m, long bitprec)
     164             : {
     165        1827 :   pari_sp av = avma;
     166        1827 :   long d = lg(Vga)-1, N, j, l, limn, prec;
     167             :   GEN LA, lj, mj, mat;
     168        1827 :   double C2 = MELLININV_CUTOFF;
     169             : 
     170        1827 :   LA = gammapoles(Vga);
     171        1827 :   N = lg(LA)-1;
     172        1827 :   lj = cgetg(N+1, t_VECSMALL);
     173        1827 :   mj = cgetg(N+1, t_VEC);
     174        4669 :   for (j = 1; j <= N; ++j)
     175             :   {
     176        2842 :     GEN L = gel(LA,j);
     177        2842 :     lj[j] = lg(L)-1;
     178        2842 :     gel(mj,j) = gsubsg(2, vecmin(L));
     179             :   }
     180        1827 :   prec = nbits2prec((long)(1+bitprec*(1+M_PI*d/C2)));
     181        1827 :   limn = ceil(2*M_LN2*bitprec/(d*dbllambertW0(C2/(M_PI*M_E))));
     182        1827 :   mat = cgetg(N+1, t_VEC);
     183        1827 :   l = limn + 2;
     184        4669 :   for (j=1; j <= N; j++)
     185             :   {
     186        2842 :     GEN C, c, mjj = gel(mj,j), pr = gen_1, t = gen_1;
     187        2842 :     long i, k, n, ljj = lj[j], lprecdl = ljj+3;
     188       12089 :     for (i=1; i <= d; i++)
     189             :     {
     190        9247 :       GEN a = gmul2n(gadd(mjj, gel(Vga,i)), -1);
     191        9247 :       GEN u = deg1pol_shallow(ghalf, a, 0);
     192        9247 :       pr = gmul(pr, ggamma(RgX_to_ser(u, lprecdl), prec));
     193        9247 :       t = gmul(t, u);
     194             :     }
     195        2842 :     c = cgetg(limn+2,t_COL);
     196        2842 :     gel(c,1) = pr;
     197      101648 :     for (n=1; n <= limn; n++)
     198       98806 :       gel(c,n+1) = gdiv(gel(c,n), RgX_translate(t, stoi(-2*n)));
     199             : 
     200        2842 :     gel(mat, j) = C = cgetg(ljj+1, t_COL);
     201        8407 :     for (k = 1; k <= ljj; k++)
     202             :     {
     203        5565 :       GEN L = cgetg(l, t_POL);
     204        5565 :       L[1] = evalsigne(1)|evalvarn(0);
     205        5565 :       for (n = 2; n < l; n++) gel(L,n) = sercoeff(gel(c,n), -k, prec);
     206        5565 :       gel(C,k) = L;
     207             :     }
     208             :     /* C[k] = \sum_n c_{j,k} t^n =: C_k(t) in Dokchitser's Algo 3.3 */
     209             :     /* Take m-th derivative of t^(-mj+2) sum_k (-ln t)^k/k! C_k(t^2) */
     210        2842 :     if (m)
     211             :     {
     212         119 :       mjj = gsubgs(mjj, 2);
     213         238 :       for (i = 1; i <= m; i++, mjj = gaddgs(mjj, 1))
     214             :       {
     215         322 :         for (k = 1; k <= ljj; k++)
     216             :         {
     217         203 :           GEN c = gel(C,k), d = RgX_shift_shallow(gmul2n(RgX_deriv(c),1), 1);
     218         203 :           c = RgX_Rg_mul(c, mjj);
     219         203 :           if (k < ljj) c = RgX_add(c, gel(C,k+1));
     220         203 :           gel(C,k) = RgX_sub(d, c);
     221             :         }
     222             :       }
     223         119 :       gel(mj,j) = gaddgs(mjj,2);
     224             :     }
     225        8407 :     for (k = 1; k <= ljj; k++)
     226             :     {
     227        5565 :       GEN c = gel(C,k);
     228        5565 :       if (k > 2) c = RgX_Rg_div(c, mpfact(k-1));
     229        5565 :       gel(C,k) = RgX_to_RgC(c, lgpol(c));
     230             :     }
     231             :   }
     232             :   /* Algo 3.3: * \phi^(m)(t) = sum_j t^m_j sum_k (-ln t)^k mat[j,k](t^2) */
     233        1827 :   return gerepilecopy(av, mkvec3(lj,RgV_neg(mj),mat));
     234             : }
     235             : 
     236             : /* Evaluate a vector considered as a polynomial using Horner. Unstable!
     237             :  * If ui != NULL, ui = 1/u, evaluate P(1/u)*u^(deg P): useful for |u|>1 */
     238             : static GEN
     239      207705 : evalvec(GEN vec, long lim, GEN u, GEN ui)
     240             : {
     241      207705 :   pari_sp ltop = avma;
     242      207705 :   GEN S = gen_0;
     243             :   long n;
     244      207705 :   lim = minss(lim, lg(vec)-1);
     245      207705 :   if (!ui)
     246       53874 :     for (n = lim; n >= 1; --n) S = gmul(u, gadd(gel(vec,n), S));
     247             :   else
     248             :   {
     249      153831 :     for (n = 1; n <= lim; ++n) S = gmul(ui, gadd(gel(vec,n), S));
     250      153831 :     S = gmul(gpowgs(u, n), S);
     251             :   }
     252      207705 :   return gerepileupto(ltop, S);
     253             : }
     254             : 
     255             : /* gammamellininvinit accessors */
     256             : static double
     257     4592582 : get_tmax(long bitprec)
     258     4592582 : { return (M_LN2 / MELLININV_CUTOFF) * bitprec ; }
     259             : static GEN
     260     9906168 : GMi_get_Vga(GEN K) { return gel(K,2); }
     261             : static long
     262     4977246 : GMi_get_m(GEN K) { return itos( gel(K,3) ); }
     263             : static GEN /* [lj,mj,mat], Kderivsmall only */
     264     5037022 : GMi_get_VS(GEN K) { return gel(K,4); }
     265             : static GEN /* [Ms,cd,A2], Kderivlarge only */
     266     9834940 : GMi_get_VL(GEN K) { return gel(K,5); }
     267             : static double
     268     4977246 : GMi_get_tmax(GEN K, long bitprec)
     269     4977246 : { return (typ(GMi_get_VS(K)) == t_INT)? -1.0 : get_tmax(bitprec); }
     270             : 
     271             : /* Compute m-th derivative of inverse Mellin at x by generalized power series
     272             :  * around x = 0; x2d = x^(2/d), x is possibly NULL (don't bother about
     273             :  * complex branches). Assume |x|^(2/d) <= tmax = M_LN2*bitprec/MELLININV_CUTOFF*/
     274             : static GEN
     275       59776 : Kderivsmall(GEN K, GEN x, GEN x2d, long bitprec)
     276             : {
     277       59776 :   pari_sp ltop = avma;
     278       59776 :   GEN Vga = GMi_get_Vga(K), VS = GMi_get_VS(K);
     279       59776 :   GEN lj = gel(VS,1), mj = gel(VS,2), mat = gel(VS,3);
     280             :   GEN d2, Lx, x2, x2i, A, S, pi;
     281       59776 :   long prec, d, N, j, k, limn, m = GMi_get_m(K);
     282             :   double Ed, xd, Wd;
     283             : 
     284       59776 :   N = lg(lj)-1; d = lg(Vga)-1;
     285       59776 :   Ed = M_LN2*bitprec / d;
     286       59776 :   xd = maxdd(M_PI*dblmodulus(x2d), 1E-13); /* pi |x|^2/d unless x tiny */
     287       59776 :   if (xd > Ed) pari_err_BUG("Kderivsmall (x2d too large)");
     288             :   /* Lemma 5.2.6 (2), a = 1 + log(Pi x^(2/d)) = log(e / xd),
     289             :    * B = log(2)*bitprec / d = Ed */
     290       59776 :   Wd = dbllambertW0( Ed / (M_E*xd) ); /* solution of w exp(w) = B exp(-a)*/
     291       59776 :   limn = (long) ceil(2*Ed/Wd);
     292       59776 :   prec = nbits2prec((long) ceil(bitprec+d*xd/M_LN2));
     293       59776 :   pi = mppi(prec);
     294       59776 :   d2 = gdivsg(d,gen_2);
     295       59776 :   if (x)
     296        1052 :     x = gmul(gtofp(x,prec), gpow(pi,d2,prec));
     297             :   else
     298       58724 :     x = gpow(gmul(gtofp(x2d,prec),pi), d2, prec);
     299             :   /* at this stage, x has been replaced by pi^(d/2) x */
     300       59776 :   x2 = gsqr(x);
     301       59776 :   Lx = gpowers(gneg(glog(x,prec)), vecsmall_max(lj));
     302       59776 :   x2i = (gcmp(gnorml2(x2), gen_1) <= 0)? NULL: ginv(x2);
     303       59776 :   S = gen_0;
     304      155957 :   for (j = 1; j <= N; ++j)
     305             :   {
     306       96181 :     long ljj = lj[j];
     307       96181 :     GEN s = gen_0;
     308      303886 :     for (k = 1; k <= ljj; k++)
     309      207705 :       s = gadd(s, gmul(gel(Lx,k), evalvec(gmael(mat,j,k), limn, x2, x2i)));
     310       96181 :     S = gadd(S, gmul(gpow(x, gel(mj,j), prec), s));
     311             :   }
     312       59776 :   A = gsubsg(m*d, sumVga(Vga));
     313       59776 :   if (!gequal0(A)) S = gmul(S, gpow(pi, gmul2n(A,-1), prec));
     314       59776 :   return gerepileupto(ltop, gtofp(S, nbits2prec(bitprec)));
     315             : }
     316             : 
     317             : /* In Klarge, we conpute K(t) as (asymptotic) * F(z), where F ~ 1 is given by
     318             :  * a continued fraction and z = Pi t^(2/d). If we take 2n terms in F (n terms
     319             :  * in Euler form), F_n(z) - F(z) is experimentally in exp(- C sqrt(n*z))
     320             :  * where C ~ 8 for d > 2 [HEURISTIC] and C = 4 (theorem) for d = 1 or d = 2
     321             :  * and vga = [0,1]. For e^(-E) absolute error, we want
     322             :  *   exp(-C sqrt(nz)) < e^-(E+a), where a ~ ln(asymptotic)
     323             :  * i.e.  2n > (E+a)^2 / t^(2/d) * 2/(C^2 Pi); C^2*Pi/2 ~ 100.5 ~ 101
     324             :  *
     325             :  * In fact, this model becomes wrong for z large: we use instead
     326             :  *
     327             :  *   exp(- sqrt(D * nz/log(z+1))) < e^-(E+a),
     328             :  * i.e.  2n > (E+a)^2 * log(1 + Pi t^(2/d))/ t^(2/d) * 2/(D Pi); */
     329             : static double
     330     4929251 : get_D(long d) { return d <= 2 ? 157. : 180.; }
     331             : /* if (abs), absolute error rather than relative */
     332             : static void
     333     4917470 : Kderivlarge_optim(GEN K, long abs, GEN t2d,GEN gcd, long *pbitprec, long *pnlim)
     334             : {
     335     4917470 :   GEN Vga = GMi_get_Vga(K), VL = GMi_get_VL(K), A2 = gel(VL,3);
     336     4917470 :   long bitprec = *pbitprec, d = lg(Vga)-1;
     337     4917470 :   const double D = get_D(d), td = dblmodulus(t2d), cd = gtodouble(gcd);
     338     4917470 :   double a, rtd, E = M_LN2*bitprec;
     339             : 
     340     4917470 :   rtd = (typ(t2d) == t_COMPLEX)? gtodouble(gel(t2d,1)): td;
     341             : 
     342             :   /* A2/2 = A, log(td) = (2/d)*log t */
     343     4917470 :   a = d*gtodouble(A2)*log2(td)/2 - (M_PI/M_LN2)*d*rtd + log2(cd); /*log2 K(t)~a*/
     344             : 
     345             :   /* if bitprec <= 0, caller should return K(t) ~ 0 */
     346     4917470 :   bitprec += 64;
     347     4917470 :   if (abs)
     348             :   {
     349     4912789 :     bitprec += ceil(a);
     350     4912789 :     if (a <= -65) E = M_LN2*bitprec; /* guarantees E <= initial E */
     351             :   }
     352     4917470 :   *pbitprec = bitprec;
     353     4917470 :   *pnlim = ceil(E*E * log2(1+M_PI*td) / (D*td));
     354     4917470 : }
     355             : 
     356             : /* Compute m-th derivative of inverse Mellin at t by continued fraction of
     357             :  * asymptotic expansion; t2d = t^(2/d). If t is NULL, "lfun" mode: don't
     358             :  * bother about complex branches + use absolute (rather than relative)
     359             :  * accuracy */
     360             : static GEN
     361     4917470 : Kderivlarge(GEN K, GEN t, GEN t2d, long bitprec0)
     362             : {
     363     4917470 :   pari_sp ltop = avma;
     364     4917470 :   GEN tdA, P, S, pi, z, Vga = GMi_get_Vga(K);
     365     4917470 :   const long d = lg(Vga)-1;
     366     4917470 :   GEN M, VL = GMi_get_VL(K), Ms = gel(VL,1), cd = gel(VL,2), A2 = gel(VL,3);
     367     4917470 :   long status, prec, nlim, m = GMi_get_m(K), bitprec = bitprec0;
     368             : 
     369     4917470 :   Kderivlarge_optim(K, !t, t2d, cd, &bitprec, &nlim);
     370     4917470 :   if (bitprec <= 0) return gen_0;
     371     4850277 :   prec = nbits2prec(bitprec);
     372     4850277 :   t2d = gtofp(t2d, prec);
     373     4850277 :   if (t)
     374        4681 :     tdA = gpow(t, gdivgs(A2,d), prec);
     375             :   else
     376     4845596 :     tdA = gpow(t2d, gdivgs(A2,2), prec);
     377     4850277 :   tdA = gmul(cd, tdA);
     378             : 
     379     4850277 :   pi = mppi(prec);
     380     4850277 :   z = gmul(pi, t2d);
     381     4850277 :   P = gmul(tdA, gexp(gmulsg(-d, z), prec));
     382     4850277 :   if (m) P = gmul(P, gpowgs(mulsr(-2, pi), m));
     383     4850277 :   M = gel(Ms,1);
     384     4850277 :   status = itos(gel(Ms,2));
     385     4850277 :   if (status == 2)
     386             :   {
     387      394303 :     if (lg(M) == 2) /* shortcut: continued fraction is constant */
     388      393974 :       S = gel(M,1);
     389             :     else
     390         329 :       S = poleval(RgV_to_RgX(M, 0), ginv(z));
     391             :   }
     392             :   else
     393             :   {
     394     4455974 :     S = contfraceval_inv(M, z, nlim/2);
     395     4455974 :     if (DEBUGLEVEL>3)
     396             :     {
     397           0 :       GEN S0 = contfraceval_inv(M, z, nlim/2 + 1);
     398           0 :       long e = gexpo(gmul(P, gabs(gsub(S,S0),0)));
     399           0 :       if (-e < bitprec0)
     400           0 :         err_printf("Kderivlarge: e = %ld, bit = %ld\n",e,bitprec0);
     401             :     }
     402     4455974 :     if (status == 1) S = gmul(S, gsubsg(1, ginv(gmul(z, pi))));
     403             :   }
     404     4850277 :   return gerepileupto(ltop, gmul(P, S));
     405             : }
     406             : 
     407             : /* Dokchitser's coefficients used for asymptotic expansion of inverse Mellin
     408             :  * 2 <= p <= min(n+1, d) */
     409             : static GEN
     410      486482 : fun_vp(long p, long n, long d, GEN SM, GEN vsinh)
     411             : {
     412      486482 :   pari_sp ltop = avma;
     413             :   long m, j, k;
     414      486482 :   GEN s = gen_0;
     415     2302206 :   for (m = 0; m <= p; ++m)
     416             :   {
     417     1815724 :     GEN pr = gen_1, s2 = gen_0, sh = gel(vsinh, d-p+1);/* (sh(x)/x)^(d-p) */
     418     1815724 :     long pm = p-m;
     419     1815724 :     for (j = m; j < p; ++j) pr = muliu(pr, d-j);
     420     4599016 :     for (k = 0; k <= pm; k+=2)
     421             :     {
     422     2783292 :       GEN e = gdiv(powuu(2*n-p+1, pm-k), mpfact(pm-k));
     423     2783292 :       s2 = gadd(s2, gmul(e, RgX_coeff(sh, k)));
     424             :     }
     425     1815724 :     s = gadd(s, gmul(gmul(gel(SM, m+1), pr), s2));
     426     1815724 :     if (gc_needed(ltop, 1)) s = gerepilecopy(ltop, s);
     427             :   }
     428      486482 :   return gerepileupto(ltop, gmul(gdivsg(-d, powuu(2*d, p)), s));
     429             : }
     430             : 
     431             : /* Asymptotic expansion of inverse Mellin, to length nlimmax. Set status = 0
     432             :  * (regular), 1 (one Hankel determinant vanishes => contfracinit will fail)
     433             :  * or 2 (same as 1, but asymptotic expansion is finite!)
     434             :  *
     435             :  * If status = 2, the asymptotic expansion is finite so return only
     436             :  * the necessary number of terms nlim <= nlimmax + d. */
     437             : static GEN
     438       11795 : Klargeinit0(GEN Vga, long nlimmax, long *status)
     439             : {
     440       11795 :   const long prec = LOWDEFAULTPREC;
     441       11795 :   const long d = lg(Vga)-1;
     442             :   long k, n, m, cnt;
     443             :   GEN pol, SM, nS1, se, vsinh, M, dk;
     444             : 
     445       11795 :   if (d==1 || (d==2 && gequal1(gabs(gsub(gel(Vga,1), gel(Vga,2)), prec))))
     446             :   { /* shortcut */
     447        9933 :     *status = 2; return mkvec(gen_1);
     448             :   }
     449             :   /* d >= 2 */
     450        1862 :   *status = 0;
     451        1862 :   pol = roots_to_pol(gneg(Vga), 0); /* deg(pol) = d */
     452        1862 :   nS1 = gpowers(gneg(RgX_coeff(pol, d-1)), d);
     453        1862 :   dk = gpowers(utoi(d), d-1);
     454        1862 :   SM = cgetg(d+3, t_VEC);
     455        9373 :   for (m = 0; m <= d; ++m)
     456             :   {
     457        7511 :     pari_sp btop = avma;
     458        7511 :     GEN s = gmul(gdivgs(gel(nS1, m+1), d), binomialuu(d, m));
     459       19978 :     for (k = 1; k <= m; ++k)
     460             :     {
     461       12467 :       GEN e = gmul(gel(nS1, m-k+1), gel(dk, k));
     462       12467 :       s = gadd(s, gmul(gmul(e, binomialuu(d-k, m-k)), RgX_coeff(pol, d-k)));
     463             :     }
     464        7511 :     gel(SM, m+1) = gerepileupto(btop, s);
     465             :   }
     466        1862 :   se = gdiv(gsinh(RgX_to_ser(pol_x(0), d+2), prec), pol_x(0));
     467        1862 :   vsinh = gpowers(se, d);
     468        1862 :   M = vectrunc_init(nlimmax + d);
     469        1862 :   vectrunc_append(M, gen_1);
     470      248110 :   for (n=2, cnt=0; (n <= nlimmax) || cnt; ++n)
     471             :   {
     472      248110 :     pari_sp btop = avma;
     473      248110 :     long p, ld = minss(d, n);
     474      248110 :     GEN s = gen_0;
     475      734592 :     for (p = 2; p <= ld; ++p)
     476      486482 :       s = gadd(s, gmul(fun_vp(p, n-1, d, SM, vsinh), gel(M, n+1-p)));
     477      248110 :     s = gerepileupto(btop, gdivgs(s, n-1));
     478      248110 :     vectrunc_append(M, s);
     479      248110 :     if (!isintzero(s))
     480             :     {
     481      248061 :       if (n >= nlimmax) break;
     482      246227 :       cnt = 0;
     483             :     }
     484             :     else
     485             :     {
     486          49 :       cnt++; *status = 1;
     487          49 :       if (cnt >= d-1) { *status = 2; setlg(M, lg(M) - (d-1)); break; }
     488             :     }
     489             :   }
     490        1862 :   return M;
     491             : }
     492             : 
     493             : /* remove trailing zeros from vector. */
     494             : static void
     495         252 : stripzeros(GEN M)
     496             : {
     497             :   long i;
     498         441 :   for(i = lg(M)-1; i >= 1; --i)
     499         441 :     if (!gequal0(gel(M, i))) break;
     500         252 :   setlg(M, i+1);
     501         252 : }
     502             : 
     503             : /* Asymptotic expansion of the m-th derivative of inverse Mellin, to length
     504             :  * nlimmax. If status = 2, the asymptotic expansion is finite so return only
     505             :  * the necessary number of terms nlim <= nlimmax + d. */
     506             : static GEN
     507       11795 : gammamellininvasymp_i(GEN Vga, long nlimmax, long m, long *status)
     508             : {
     509       11795 :   pari_sp ltop = avma;
     510             :   GEN M, A, Aadd;
     511             :   long d, i, nlim, n;
     512             : 
     513       11795 :   M = Klargeinit0(Vga, nlimmax, status);
     514       11795 :   if (!m) return gerepilecopy(ltop, M);
     515         252 :   d = lg(Vga)-1;
     516             :   /* half the exponent of t in asymptotic expansion. */
     517         252 :   A = gdivgs(gaddsg(1-d, sumVga(Vga)), 2*d);
     518         252 :   if (*status == 2) M = shallowconcat(M, zerovec(m));
     519         252 :   nlim = lg(M)-1;
     520         252 :   Aadd = gdivgs(stoi(2-d), 2*d); /* (1/d) - (1/2) */
     521         532 :   for (i = 1; i <= m; i++, A = gadd(A,Aadd))
     522        8022 :     for (n = nlim-1; n >= 1; --n)
     523       15484 :       gel(M, n+1) = gsub(gel(M, n+1),
     524        7742 :                          gmul(gel(M, n), gsub(A, gdivgs(stoi(n-1), d))));
     525         252 :   stripzeros(M);
     526         252 :   return gerepilecopy(ltop, M);
     527             : }
     528             : GEN
     529          14 : gammamellininvasymp(GEN Vga, long nlimmax, long m)
     530             : { long status;
     531          14 :   if (!is_vec_t(typ(Vga))) pari_err_TYPE("gammamellininvinit",Vga);
     532          14 :   return gammamellininvasymp_i(Vga, nlimmax, m, &status);
     533             : }
     534             : 
     535             : /* Does the continued fraction of the asymptotic expansion M at oo of inverse
     536             :  * Mellin transform attached to Vga have zero Hankel determinants ? */
     537             : static long
     538        1820 : ishankelspec(GEN Vga, GEN M)
     539             : {
     540        1820 :   long status, i, d = lg(Vga)-1;
     541             : 
     542        1820 :   if (d == 5 || d == 7)
     543             :   { /* known bad cases: a x 5 and a x 7 */
     544          70 :     GEN v1 = gel(Vga, 1);
     545         294 :     for (i = 2; i <= d; ++i)
     546         238 :       if (!gequal(gel(Vga,i), v1)) break;
     547          70 :     if (i > d) return 1;
     548             :   }
     549        1764 :   status = 0;
     550             :   /* Heuristic: if 6 first terms in contfracinit don't fail, assume it's OK */
     551        1764 :   pari_CATCH(e_INV) { status = 1; }
     552        1764 :   pari_TRY { contfracinit(M, minss(lg(M)-2,6)); }
     553        1764 :   pari_ENDCATCH; return status;
     554             : }
     555             : 
     556             : /* Initialize data for computing m-th derivative of inverse Mellin */
     557             : GEN
     558       11781 : gammamellininvinit(GEN Vga, long m, long bitprec)
     559             : {
     560       11781 :   pari_sp ltop = avma;
     561             :   GEN A2, M, VS, VL, cd;
     562       11781 :   long d = lg(Vga)-1, status;
     563       11781 :   const double C2 = MELLININV_CUTOFF, D = get_D(d);
     564       11781 :   double E = M_LN2*bitprec, tmax = get_tmax(bitprec); /* = E/C2 */
     565       11781 :   const long nlimmax = ceil(E*log2(1+M_PI*tmax)*C2/D);
     566             : 
     567       11781 :   if (!is_vec_t(typ(Vga))) pari_err_TYPE("gammamellininvinit",Vga);
     568       11781 :   A2 = gaddsg(m*(2-d) + 1-d, sumVga(Vga));
     569       11781 :   cd = (d <= 2)? gen_2: gsqrt(gdivgs(int2n(d+1), d), nbits2prec(bitprec));
     570             :   /* if in Klarge, we have |t| > tmax = E/C2, thus nlim < E*C2/D. */
     571       11781 :   M = gammamellininvasymp_i(Vga, nlimmax, m, &status);
     572       11781 :   if (status == 2)
     573             :   {
     574        9954 :     tmax = -1.; /* only use Klarge */
     575        9954 :     VS = gen_0;
     576             :   }
     577             :   else
     578             :   {
     579        1827 :     long prec = nbits2prec((4*bitprec)/3);
     580        1827 :     VS = Kderivsmallinit(Vga, m, bitprec);
     581        1827 :     if (status == 0 && ishankelspec(Vga, M)) status = 1;
     582        1827 :     if (status == 1)
     583             :     { /* a Hankel determinant vanishes => contfracinit is undefined.
     584             :          So compute K(t) / (1 - 1/(pi^2*t)) instead of K(t)*/
     585          63 :       GEN t = ginv(mppi(prec));
     586             :       long i;
     587        9713 :       for (i = 2; i < lg(M); ++i)
     588        9650 :         gel(M, i) = gadd(gel(M, i), gmul(gel(M, i-1), t));
     589             :     }
     590             :     else
     591        1764 :       M = RgC_gtofp(M, prec); /* convert from rationals to t_REAL: faster */
     592        1827 :     M = contfracinit(M, lg(M)-2);
     593             :   }
     594       11781 :   VL = mkvec3(mkvec2(M, stoi(status)), cd, A2);
     595       11781 :   return gerepilecopy(ltop, mkvec5(dbltor(tmax), Vga, stoi(m), VS, VL));
     596             : }
     597             : 
     598             : /* Compute m-th derivative of inverse Mellin at s2d = s^(d/2) using
     599             :  * initialization data. Use Taylor expansion at 0 for |s2d| < tmax, and
     600             :  * asymptotic expansion at oo otherwise. WARNING: assume that accuracy
     601             :  * has been increased according to tmax by the CALLING program. */
     602             : GEN
     603     4971513 : gammamellininvrt(GEN K, GEN s2d, long bitprec)
     604             : {
     605     4971513 :   if (dblmodulus(s2d) < GMi_get_tmax(K, bitprec))
     606       58724 :     return Kderivsmall(K, NULL, s2d, bitprec);
     607             :   else
     608     4912789 :     return Kderivlarge(K, NULL, s2d, bitprec);
     609             : }
     610             : 
     611             : /* Compute inverse Mellin at s. K from gammamellininv OR a Vga, in which
     612             :  * case the initialization data is computed. */
     613             : GEN
     614        5733 : gammamellininv(GEN K, GEN s, long m, long bitprec)
     615             : {
     616        5733 :   pari_sp av = avma;
     617             :   GEN z, s2d;
     618             :   long d;
     619        5733 :   if (!is_vec_t(typ(K))) pari_err_TYPE("gammamellininv",K);
     620        5733 :   if (lg(K) != 6 || !is_vec_t(typ(GMi_get_Vga(K))))
     621          14 :     K = gammamellininvinit(K, m, bitprec);
     622        5733 :   d = lg(GMi_get_Vga(K))-1;
     623        5733 :   s2d = gpow(s, gdivgs(gen_2, d), nbits2prec(bitprec));
     624        5733 :   if (dblmodulus(s2d) < GMi_get_tmax(K, bitprec))
     625        1052 :     z = Kderivsmall(K, s, s2d, bitprec);
     626             :   else
     627        4681 :     z = Kderivlarge(K, s, s2d, bitprec);
     628        5733 :   return gerepileupto(av, z);
     629             : }

Generated by: LCOV version 1.13