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

Generated by: LCOV version 1.11