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

Generated by: LCOV version 1.13