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

Generated by: LCOV version 1.11