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.16.2 lcov report (development 29115-f22e516b23) Lines: 349 354 98.6 %
Date: 2024-04-22 08:08:03 Functions: 40 40 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; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : #include "pari.h"
      16             : #include "paripriv.h"
      17             : 
      18             : #define DEBUGLEVEL DEBUGLEVEL_gammamellininv
      19             : 
      20             : /*******************************************************************/
      21             : /*               Computation of inverse Mellin                     */
      22             : /*               transforms of gamma products.                     */
      23             : /*******************************************************************/
      24             : /* Handle complex Vga whose sum is real */
      25             : static GEN
      26       17874 : sumVga(GEN Vga) { return real_i(vecsum(Vga)); }
      27             : 
      28             : /* ac != 0 */
      29             : static double
      30      562938 : lemma526_i(double ac, double c, double t, double B)
      31             : {
      32      562938 :   double D = -B/ac; /* sgn(t) = sgn(a) = - sgn(D) */
      33      562938 :   if (D <= 0)
      34             :   {
      35      489324 :     if (D > -100)
      36             :     {
      37       67291 :       D = -exp(D) / t;
      38       67291 :       if (D < - 1/M_E) return 0;
      39       66962 :       D = dbllambertW_1(D);
      40             :     }
      41             :     else
      42             :     { /* avoid underflow, use asymptotic expansion */
      43      422033 :       double U = D - log(t);
      44      422033 :       D = U - log(-U);
      45             :     }
      46      488995 :     return pow(maxdd(t, -t * D), c);
      47             :   }
      48             :   else
      49             :   {
      50       73614 :     if (D < 100)
      51        9702 :       D = dbllambertW0(-exp(D) / t);
      52             :     else
      53             :     { /* avoid overflow, use asymptotic expansion */
      54       63912 :       double U = D - log(-t);
      55       63912 :       D = U - log(U);
      56             :     }
      57       73614 :     return pow(-t * D, c);
      58             :   }
      59             : }
      60             : /* b > 0, c > 0; solve x^a exp(-b x^(1/c)) < e^(-B) for x >= 0 */
      61             : double
      62          14 : dbllemma526(double a, double b, double c, double B)
      63             : {
      64             :   double ac;
      65          14 :   if (!a) return B <= 0? 0: pow(B/b, c);
      66          14 :   ac = a*c; if (B < 0) B = 1e-9;
      67          14 :   return lemma526_i(ac, c, ac/b, B);
      68             : }
      69             : /* Same, special case b/c = 2Pi, the only one needed: for c = d/2 */
      70             : double
      71     2354096 : dblcoro526(double a, double c, double B)
      72             : {
      73     2354096 :   if (!a) return B <= 0? 0: pow(B/(2*M_PI*c), c);
      74      562924 :   if (B < 0) B = 1e-9;
      75      562924 :   return lemma526_i(a*c, c, a/(2*M_PI), B);
      76             : }
      77             : 
      78             : static const double MELLININV_CUTOFF = 121.; /* C*C */
      79             : 
      80             : /* x real */
      81             : static GEN
      82        9457 : RMOD2(GEN x) { return gsub(x, gmul2n(gdiventgs(x,2), 1)); }
      83             : /* x real or complex, return canonical representative for x mod 2Z */
      84             : static GEN
      85        9457 : MOD2(GEN x)
      86        9457 : { return typ(x) == t_COMPLEX? mkcomplex(RMOD2(gel(x,1)), gel(x,2)): RMOD2(x); }
      87             : static GEN
      88        2898 : RgV_MOD2(GEN x)
      89       12355 : { pari_APPLY_same(MOD2(gel(x,i))); }
      90             : 
      91             : /* classes of poles of the gamma factor mod 2Z, sorted by increasing
      92             :  * Re(s) mod 2 (in [0,2[).*/
      93             : static GEN
      94        2898 : gammapoles(GEN Vga, long *pdV, long bit)
      95             : {
      96        2898 :   long i, m, emax, l = lg(Vga);
      97        2898 :   GEN P, B = RgV_MOD2(Vga), V = cgetg(l, t_VEC);
      98        2898 :   P = gen_indexsort(B, (void*)lexcmp, cmp_nodata);
      99        7994 :   for (i = m = 1; i < l;)
     100             :   {
     101        5096 :     GEN u = gel(B, P[i]);
     102             :     long k;
     103        9457 :     for(k = i+1; k < l; k++)
     104             :     {
     105        6559 :       GEN v = gsub(u, gel(B, P[k]));
     106        6559 :       if (!gequal0(v) && (!isinexactreal(v) || gexpo(v) > -bit)) break;
     107             :     }
     108        5096 :     gel(V, m++) = vecslice(P,i,k-1);
     109        5096 :     i = k;
     110             :   }
     111        2898 :   setlg(V, m); emax = 0;
     112        7994 :   for (i = 1; i < m; i++)
     113             :   {
     114        5096 :     long j, e = 0, li = lg(gel(V,i))-1;
     115        5096 :     GEN b = gel(B, gel(V,i)[1]);
     116       15526 :     for (j = 1; j < m; j++)
     117       10430 :       if (j != i) e -= gexpo(gsub(gel(B, gel(V,j)[1]), b));
     118        5096 :     emax = maxss(emax, e * li);
     119             :   }
     120        7994 :   for (i = 1; i < m; i++) gel(V,i) = vecpermute(Vga, gel(V,i));
     121        2898 :   *pdV = emax; return V;
     122             : }
     123             : 
     124             : static GEN
     125      474390 : sercoeff(GEN x, long n, long prec)
     126             : {
     127      474390 :   long N = n - valser(x);
     128      474390 :   return (N < 0)? gen_0: gprec_wtrunc(gel(x, N+2), prec);
     129             : }
     130             : 
     131             : /* prod_i Gamma(s/2 + (m+LA[i])/2), set t *= prod_i (s/2 + (m+LA[i])/2) */
     132             : static GEN
     133       10430 : get_gamma(GEN *pt, GEN LA, GEN m, int round, long precdl, long prec)
     134             : {
     135       10430 :   long i, l = lg(LA);
     136       10430 :   GEN pr = NULL, t = *pt;
     137       27825 :   for (i = 1; i < l; i++)
     138             :   {
     139       17395 :     GEN u, g, a = gmul2n(gadd(m, gel(LA,i)), -1);
     140       17395 :     if (round) a = ground(a);
     141       17395 :     u = deg1pol_shallow(ghalf, a, 0);
     142       17395 :     g = ggamma(RgX_to_ser(u, precdl), prec);
     143       17395 :     pr = pr? gmul(pr, g): g;
     144       17395 :     t = t? gmul(t, u): u;
     145             :   }
     146       10430 :   *pt = t; return pr;
     147             : }
     148             : /* generalized power series expansion of inverse Mellin around x = 0;
     149             :  * m-th derivative */
     150             : static GEN
     151        2898 : Kderivsmallinit(GEN ldata, GEN Vga, long m, long bit)
     152             : {
     153        2898 :   const double C2 = MELLININV_CUTOFF;
     154        2898 :   long prec2, N, j, l, dLA, limn, d = lg(Vga)-1;
     155             :   GEN piA, LA, L, M, mat;
     156             : 
     157        2898 :   LA = gammapoles(Vga, &dLA, bit); N = lg(LA)-1;
     158        2898 :   prec2 = nbits2prec(dLA + bit * (1 + M_PI*d/C2));
     159             : #if BITS_IN_LONG == 32
     160         414 :   prec2 += odd(prec2lg(prec2)) ? EXTRAPRECWORD: 0;
     161             : #endif
     162        2898 :   if (ldata) Vga = ldata_get_gammavec(ldata_newprec(ldata, prec2));
     163        2898 :   L = cgetg(N+1, t_VECSMALL);
     164        2898 :   M = cgetg(N+1, t_VEC);
     165        2898 :   mat = cgetg(N+1, t_VEC);
     166        2898 :   limn = ceil(2*M_LN2*bit / (d * dbllambertW0(C2/(M_PI*M_E*d))));
     167        2898 :   l = limn + 2;
     168        7994 :   for (j = 1; j <= N; j++)
     169             :   {
     170        5096 :     GEN S = gel(LA,j);
     171        5096 :     GEN C, c, mj, G = NULL, t = NULL, tj = NULL;
     172        5096 :     long i, k, n, jj, lj = L[j] = lg(S)-1, precdl = lj+3;
     173             : 
     174        5096 :     gel(M,j) = mj = gsubsg(2, gel(S, vecindexmin(real_i(S))));
     175       15526 :     for (jj = 1; jj <= N; jj++)
     176             :     {
     177             :       GEN g;
     178       10430 :       if (jj == j) /* poles come from this class only */
     179        5096 :         g = get_gamma(&tj, gel(LA,jj), mj, 1, precdl, prec2);
     180             :       else
     181        5334 :         g = get_gamma(&t, gel(LA,jj), mj, 0, precdl, prec2);
     182       10430 :       G = G? gmul(G, g): g;
     183             :     }
     184        5096 :     c = cgetg(limn+2,t_COL); gel(c,1) = G;
     185      271712 :     for (n=1; n <= limn; n++)
     186             :     {
     187      266616 :       GEN A = utoineg(2*n), T = RgX_translate(tj, A);
     188             :       /* T = exact polynomial, may vanish at 0 (=> pole in c[n+1]) */
     189      266616 :       if (t) T = RgX_mul(T, RgX_translate(t, A)); /* no pole here */
     190      266616 :       gel(c,n+1) = gdiv(gel(c,n), T);
     191             :     }
     192        5096 :     gel(mat, j) = C = cgetg(lj+1, t_COL);
     193       14553 :     for (k = 1; k <= lj; k++)
     194             :     {
     195        9457 :       GEN L = cgetg(l, t_POL);
     196      483847 :       for (n = 2; n < l; n++) gel(L,n) = sercoeff(gel(c,n), -k, prec2);
     197        9457 :       L[1] = evalsigne(1)|evalvarn(0); gel(C,k) = L;
     198             :     }
     199             :     /* C[k] = \sum_n c_{j,k} t^n =: C_k(t) in Dokchitser's Algo 3.3
     200             :      * m-th derivative of t^(-M+2) sum_k (-ln t)^k/k! C_k(t^2) */
     201        5096 :     if (m)
     202             :     {
     203         119 :       mj = gsubgs(mj, 2);
     204         238 :       for (i = 1; i <= m; i++, mj = gaddgs(mj,1))
     205         322 :         for (k = 1; k <= lj; k++)
     206             :         {
     207         203 :           GEN c = gel(C,k), d = RgX_shift_shallow(gmul2n(RgX_deriv(c),1), 1);
     208         203 :           c = RgX_Rg_mul(c, mj);
     209         203 :           if (k < lj) c = RgX_add(c, gel(C,k+1));
     210         203 :           gel(C,k) = RgX_sub(d, c);
     211             :         }
     212         119 :       gel(M,j) = gaddgs(mj,2);
     213             :     }
     214       14553 :     for (k = 1; k <= lj; k++)
     215             :     {
     216        9457 :       GEN c = gel(C,k);
     217        9457 :       if (k > 2) c = RgX_Rg_div(c, mpfact(k-1));
     218        9457 :       gel(C,k) = RgX_to_RgC(c, lgpol(c));
     219             :     }
     220             :   }
     221             :   /* Algo 3.3: * \phi^(m)(t) = sum_j t^m_j sum_k (-ln t)^k mat[j,k](t^2) */
     222        2898 :   piA = gsubsg(m*d, sumVga(Vga));
     223        2898 :   if (!gequal0(piA)) piA = powPis(gmul2n(piA,-1), prec2);
     224        2898 :   return mkvec5(L, RgV_neg(M), mat, mkvecsmall(prec2), piA);
     225             : }
     226             : 
     227             : /* Evaluate a vector considered as a polynomial using Horner. */
     228             : static GEN
     229      257425 : evalvec(GEN vec, long N, GEN u)
     230             : {
     231      257425 :   GEN S = gen_0;
     232             :   long n;
     233      257425 :   N = minss(N, lg(vec)-1);
     234     8940796 :   for (n = N; n >= 1; n--) S = gmul(u, gadd(gel(vec,n), S));
     235      257427 :   return S;
     236             : }
     237             : 
     238             : /* gammamellininvinit accessors */
     239             : static double
     240     4457006 : get_tmax(long bitprec)
     241     4457006 : { return (M_LN2 / MELLININV_CUTOFF) * bitprec ; }
     242             : static GEN
     243        5726 : GMi_get_Vga(GEN K) { return gel(K,2); }
     244             : static long
     245    11206795 : GMi_get_degree(GEN K) { return lg(gel(K,2))-1; }
     246             : static long
     247     5564259 : GMi_get_m(GEN K) { return itos( gel(K,3) ); }
     248             : static GEN /* [lj,mj,mat,prec2], Kderivsmall only */
     249     5711997 : GMi_get_VS(GEN K) { return gel(K,4); }
     250             : /* K[5] = [Ms,cd,A2], Kderivlarge only */
     251             : static long/*Kderivlarge*/
     252     5564124 : GMi_get_status(GEN K) { return itos(gmael3(K,5,1,2)); }
     253             : static GEN/*Kderivlarge*/
     254     5564101 : GMi_get_M(GEN K) { return gmael3(K,5,1,1); }
     255             : static GEN/*Kderivlarge*/
     256     5564100 : GMi_get_cd(GEN K) { return gmael(K,5,2); }
     257             : static GEN/*Kderivlarge*/
     258    11126943 : GMi_get_A2(GEN K) { return gmael(K,5,3); }
     259             : 
     260             : static double
     261     5638175 : GMi_get_tmax(GEN K, long bitprec)
     262     5638175 : { return (typ(GMi_get_VS(K)) == t_INT)? -1.0 : get_tmax(bitprec); }
     263             : 
     264             : /* Compute m-th derivative of inverse Mellin at x by generalized power series
     265             :  * around x = 0; x2d = x^(2/d), x is possibly NULL (don't bother about
     266             :  * complex branches). Assume |x|^(2/d) <= tmax = M_LN2*bitprec/MELLININV_CUTOFF*/
     267             : static GEN
     268       73971 : Kderivsmall(GEN K, GEN x, GEN x2d, long bitprec)
     269             : {
     270       73971 :   GEN VS = GMi_get_VS(K), L = gel(VS,1), M = gel(VS,2), mat = gel(VS,3);
     271       73971 :   GEN d2, Lx, x2, S, pi, piA = gel(VS,5);
     272       73971 :   long j, k, prec = gel(VS,4)[1], d = GMi_get_degree(K), limn, N = lg(L)-1;
     273       73971 :   double xd, Wd, Ed = M_LN2*bitprec / d;
     274             : 
     275       73971 :   xd = maxdd(M_PI*dblmodulus(x2d), 1E-13); /* pi |x|^2/d unless x tiny */
     276       73971 :   if (xd > Ed) pari_err_BUG("Kderivsmall (x2d too large)");
     277             :   /* Lemma 5.2.6 (2), a = 1 + log(Pi x^(2/d)) = log(e / xd),
     278             :    * B = log(2)*bitprec / d = Ed */
     279       73971 :   Wd = dbllambertW0( Ed / (M_E*xd) ); /* solution of w exp(w) = B exp(-a)*/
     280       73971 :   limn = (long) ceil(2*Ed/Wd);
     281       73971 :   pi = mppi(prec);
     282       73971 :   d2 = gdivsg(d,gen_2);
     283       73971 :   if (x)
     284        1122 :     x = gmul(gtofp(x,prec), gpow(pi,d2,prec));
     285             :   else
     286       72849 :     x = gpow(gmul(gtofp(x2d,prec),pi), d2, prec);
     287             :   /* at this stage, x has been replaced by pi^(d/2) x */
     288       73971 :   x2 = gsqr(x);
     289       73972 :   Lx = gpowers(gneg(glog(x,prec)), vecsmall_max(L));
     290       73971 :   S = gen_0;
     291      201398 :   for (j = 1; j <= N; ++j)
     292             :   {
     293      127426 :     long lj = L[j];
     294      127426 :     GEN s = gen_0;
     295      384852 :     for (k = 1; k <= lj; k++)
     296      257425 :       s = gadd(s, gmul(gel(Lx,k), evalvec(gmael(mat,j,k), limn, x2)));
     297      127427 :     S = gadd(S, gmul(gpow(x, gel(M,j), prec), s));
     298             :   }
     299       73972 :   if (!gequal0(piA)) S = gmul(S, piA);
     300       73971 :   return S;
     301             : }
     302             : 
     303             : /* In Klarge, we conpute K(t) as (asymptotic) * F(z), where F ~ 1 is given by
     304             :  * a continued fraction and z = Pi t^(2/d). If we take 2n terms in F (n terms
     305             :  * in Euler form), F_n(z) - F(z) is experimentally in exp(- C sqrt(n*z))
     306             :  * where C ~ 8 for d > 2 [HEURISTIC] and C = 4 (theorem) for d = 1 or d = 2
     307             :  * and vga = [0,1]. For e^(-E) absolute error, we want
     308             :  *   exp(-C sqrt(nz)) < e^-(E+a), where a ~ ln(asymptotic)
     309             :  * i.e.  2n > (E+a)^2 / t^(2/d) * 2/(C^2 Pi); C^2*Pi/2 ~ 100.5 ~ 101
     310             :  *
     311             :  * In fact, this model becomes wrong for z large: we use instead
     312             :  *
     313             :  *   exp(- sqrt(D * nz/log(z+1))) < e^-(E+a),
     314             :  * i.e.  2n > (E+a)^2 * log(1 + Pi t^(2/d))/ t^(2/d) * 2/(D Pi); */
     315             : static double
     316     5578712 : get_D(long d) { return d <= 2 ? 157. : 180.; }
     317             : /* if (abs), absolute error rather than relative */
     318             : static void
     319     5563969 : Kderivlarge_optim(GEN K, int abs, GEN t2d, double cd, long *pbitprec, long *pnlim)
     320             : {
     321     5563969 :   GEN A2 = GMi_get_A2(K);
     322     5563911 :   long bitprec = *pbitprec, d = GMi_get_degree(K);
     323     5563979 :   const double D = get_D(d), td = dblmodulus(t2d);
     324     5563905 :   double a, rtd, E = M_LN2*bitprec;
     325             : 
     326             :   /* t = 0 can happen with finite continued fraction or easyvga */
     327     5563905 :   if (!td) { *pnlim = 0; return; }
     328     5563898 :   rtd = (typ(t2d) == t_COMPLEX)? gtodouble(gel(t2d,1)): td;
     329             :   /* A2/2 = A, log(td) = (2/d)*log t */
     330     5563898 :   a = d*gtodouble(A2)*log2(td)/2 - (M_PI/M_LN2)*d*rtd + log2(cd);/*log2 K(t)~a*/
     331             :   /* if bitprec <= 0, caller should return K(t) ~ 0 */
     332     5563530 :   bitprec += 64;
     333     5563530 :   if (abs)
     334             :   {
     335     5558856 :     bitprec += ceil(a);
     336     5558856 :     if (a <= -65) E = M_LN2*bitprec; /* guarantees E <= initial E */
     337             :   }
     338     5563530 :   *pbitprec = bitprec;
     339     5563530 :   *pnlim = ceil(E*E * log2(1+M_PI*td) / (D*td));
     340             : }
     341             : 
     342             : /* Compute m-th derivative of inverse Mellin at t by continued fraction of
     343             :  * asymptotic expansion; t2d = t^(2/d). If t is NULL, "lfun" mode: don't
     344             :  * bother about complex branches + use absolute (rather than relative)
     345             :  * accuracy */
     346             : static GEN
     347     5564201 : Kderivlarge(GEN K, GEN t, GEN t2d, long bitprec0)
     348             : {
     349             :   GEN tdA, P, S, pi, z;
     350     5564201 :   const long d = GMi_get_degree(K);
     351     5564126 :   GEN M = GMi_get_M(K), cd = GMi_get_cd(K), A2 = GMi_get_A2(K);
     352     5564124 :   long prec, nlim, status = GMi_get_status(K), m = GMi_get_m(K), bitprec = bitprec0;
     353             : 
     354     5564263 :   Kderivlarge_optim(K, !t, t2d, gtodouble(cd), &bitprec, &nlim);
     355     5563680 :   if (bitprec <= 0) return gen_0;
     356     5478546 :   prec = nbits2prec(bitprec);
     357     5478557 :   t2d = gtofp(t2d, prec);
     358     5478878 :   if (t)
     359        4702 :     tdA = gpow(t, gdivgu(A2,d), prec);
     360             :   else
     361     5474176 :     tdA = gpow(t2d, gdivgu(A2,2), prec);
     362     5478763 :   pi = mppi(prec); z = gmul(pi, t2d);
     363     5478947 :   P = gmul(gmul(cd, tdA), gexp(gmulsg(-d, z), prec));
     364     5479095 :   if (m) P = gmul(P, gpowgs(mulsr(-2, pi), m));
     365     5479074 :   if (status == 2) /* finite continued fraction */
     366     1193520 :     S = (lg(M) == 2)? gel(M,1): poleval(M, ginv(z));
     367             :   else
     368             :   {
     369     4285554 :     S = contfraceval_inv(M, z, nlim/2);
     370     4285658 :     if (DEBUGLEVEL>3)
     371             :     {
     372           0 :       GEN S0 = contfraceval_inv(M, z, minss(nlim/2 + 1, minss(lg(gel(M, 1)) - 1, lg(gel(M, 2)))));
     373           0 :       long e = gexpo(gmul(P, gsub(S,S0)));
     374           0 :       if (-e < bitprec0)
     375           0 :         err_printf("Kderivlarge: e = %ld, bit = %ld\n",e,bitprec0);
     376             :     }
     377     4285658 :     if (status == 1) S = gmul(S, gsubsg(1, ginv(gmul(z, pi))));
     378             :   }
     379     5479177 :   return gmul(P, S);
     380             : }
     381             : 
     382             : /* Dokchitser's coefficients used for asymptotic expansion of inverse Mellin
     383             :  * 2 <= p <= min(n+1, d), c = 2n-p+1; sh = (sh(x)/x)^(d-p) */
     384             : static GEN
     385      914891 : vp(long p, long c, GEN SMd, GEN sh)
     386             : {
     387      914891 :   GEN s, ve = cgetg(p+2, t_VEC);
     388             :   long m, j, k;
     389             : 
     390      914891 :   gel(ve,1) = gen_1; gel(ve,2) = utoipos(c);
     391     2678700 :   for (j = 2; j <= p; j++) gel(ve,j+1) = gdivgs(gmulgu(gel(ve,j), c), j);
     392      914891 :   s = gel(SMd, 1);
     393     3593591 :   for (m = 1; m <= p; m++)
     394             :   {
     395     2678700 :     GEN t, c = gel(SMd, m+1);
     396     2678700 :     if (gequal0(c)) continue;
     397     1402882 :     t = gel(ve, m+1);
     398     2813016 :     for (k = 1; k <= m/2; k++)
     399     1410134 :       t = gadd(t, gmul(gel(ve, m-2*k+1), RgX_coeff(sh, k)));
     400     1402882 :     s = gadd(s, gmul(c, t));
     401             :   }
     402      914891 :   return s;
     403             : }
     404             : 
     405             : static GEN
     406        5705 : get_SM(GEN Vga)
     407             : {
     408        5705 :   long k, m, d = lg(Vga)-1;
     409        5705 :   GEN pol, nS1, SM, C, t = vecsum(Vga);
     410             : 
     411        5705 :   pol = roots_to_pol(gmulgs(Vga, -d), 0); /* deg(pol) = d */
     412        5705 :   SM = cgetg(d+2, t_VEC); gel(SM,1) = gen_1;
     413        5705 :   if (gequal0(t))
     414             :   { /* shortcut */
     415        9688 :     for (m = 1; m <= d; m++) gel(SM,m+1) = gel(pol,d+2-m);
     416        2436 :     return SM;
     417             :   }
     418        3269 :   nS1 = gpowers(gneg(t), d); C = matpascal(d);
     419       14336 :   for (m = 1; m <= d; m++)
     420             :   {
     421       11067 :     pari_sp av = avma;
     422       11067 :     GEN s = gmul(gel(nS1, m+1), gcoeff(C, d+1, m+1));
     423       37338 :     for (k = 1; k <= m; k++)
     424             :     {
     425       26271 :       GEN e = gmul(gel(nS1, m-k+1), gcoeff(C, d-k+1, m-k+1));
     426       26271 :       s = gadd(s, gmul(e, RgX_coeff(pol, d-k)));
     427             :     }
     428       11067 :     gel(SM, m+1) = gerepileupto(av, s);
     429             :   }
     430        3269 :   return SM;
     431             : }
     432             : 
     433             : static GEN
     434        5705 : get_SMd(GEN Vga)
     435             : {
     436        5705 :   GEN M, SM = get_SM(Vga);
     437        5705 :   long p, m, d = lg(Vga)-1;
     438             : 
     439        5705 :   M = cgetg(d, t_VEC);
     440       18319 :   for (p = 2; p <= d; p++)
     441             :   {
     442       12614 :     GEN a = gen_1, c;
     443       12614 :     long D = d - p;
     444       12614 :     gel(M, p-1) = c = cgetg(p+2, t_COL);
     445       12614 :     gel(c, 1) = gel(SM, p+1);
     446       49238 :     for (m = 1; m <= p; m++)
     447             :     {
     448       36624 :       a = muliu(a, D + m);
     449       36624 :       gel(c, m+1) = gmul(gel(SM, p-m+1), a);
     450             :     }
     451             :   }
     452        5705 :   return M;
     453             : }
     454             : 
     455             : /* Asymptotic expansion of inverse Mellin, to length nlimmax. Set status = 0
     456             :  * (regular), 1 (one Hankel determinant vanishes => contfracinit will fail)
     457             :  * or 2 (same as 1, but asymptotic expansion is finite!)
     458             :  *
     459             :  * If status = 2, the asymptotic expansion is finite so return only
     460             :  * the necessary number of terms nlim <= nlimmax + d. */
     461             : static GEN
     462       17517 : Klargeinit(GEN Vga, long nlimmax, long *status, long prec)
     463             : {
     464       17517 :   long d = lg(Vga) - 1, p, n, cnt;
     465             :   GEN M, SMd, se, vsinh, vd;
     466             : 
     467       17517 :   if (Vgaeasytheta(Vga)) { *status = 2; return mkvec(gen_1); }
     468             :   /* d >= 2 */
     469        5705 :   *status = 0;
     470        5705 :   if (prec) prec += nbits2extraprec((prec >> 1) + BITS_IN_LONG);
     471        5705 :   SMd = get_SMd(Vga);
     472        5705 :   se = gsinh(RgX_to_ser(pol_x(0), d+2), 0); setvalser(se,0);
     473        5705 :   se = gdeflate(se, 0, 2); /* se(x^2) = sinh(x)/x */
     474        5705 :   vsinh = gpowers(se, d);
     475        5705 :   vd = gpowers(utoipos(2*d), d);
     476        5705 :   M = cgetg(nlimmax + d + 1, t_VEC); gel(M,1) = gen_1;
     477      411626 :   for (n = 2, cnt = 0; n <= nlimmax || cnt; n++)
     478             :   {
     479      411626 :     pari_sp av = avma;
     480      411626 :     long ld = minss(d, n);
     481      411626 :     GEN s = gen_0;
     482     1326517 :     for (p = 2; p <= ld; p++)
     483             :     {
     484      914891 :       GEN z = vp(p, 2*n-1-p, gel(SMd, p-1), gel(vsinh, d-p+1));
     485      914891 :       s = gadd(s, gmul(gdiv(z, gel(vd, p+1)), gel(M, n+1-p)));
     486             :     }
     487      411626 :     if (prec && !isinexact(s)) s = gtofp(s, prec);
     488      411626 :     gel(M,n) = s = gerepileupto(av, gdivgs(s, 1-n));
     489      411626 :     if (gequal0(s))
     490             :     {
     491          56 :       cnt++; *status = 1;
     492          56 :       if (cnt >= d-1) { *status = 2; n -= d-2; break; }
     493             :     }
     494             :     else
     495             :     {
     496      411570 :       if (n >= nlimmax) { n++; break; }
     497      405893 :       cnt = 0;
     498             :     }
     499             :   }
     500        5705 :   setlg(M, n); return M;
     501             : }
     502             : 
     503             : /* remove trailing zeros from vector. */
     504             : static void
     505         252 : stripzeros(GEN M)
     506             : {
     507             :   long i;
     508         441 :   for(i = lg(M)-1; i >= 1; --i)
     509         441 :     if (!gequal0(gel(M, i))) break;
     510         252 :   setlg(M, i+1);
     511         252 : }
     512             : 
     513             : /* Asymptotic expansion of the m-th derivative of inverse Mellin, to length
     514             :  * nlimmax. If status = 2, the asymptotic expansion is finite so return only
     515             :  * the necessary number of terms nlim <= nlimmax + d. */
     516             : static GEN
     517       14745 : gammamellininvasymp_i(GEN Vga, long nlimmax, long m, long *status, long prec)
     518             : {
     519             :   GEN M, A, Aadd;
     520             :   long d, i, nlim, n;
     521             : 
     522       14745 :   M = Klargeinit(Vga, nlimmax, status, prec);
     523       14745 :   if (!m) return M;
     524         252 :   d = lg(Vga)-1;
     525             :   /* half the exponent of t in asymptotic expansion. */
     526         252 :   A = gdivgu(gaddsg(1-d, sumVga(Vga)), 2*d);
     527         252 :   if (*status == 2) M = shallowconcat(M, zerovec(m));
     528         252 :   nlim = lg(M)-1;
     529         252 :   Aadd = sstoQ(2-d, 2*d); /* (1/d) - (1/2) */
     530         532 :   for (i = 1; i <= m; i++, A = gadd(A,Aadd))
     531        8022 :     for (n = nlim-1; n >= 1; --n)
     532       15484 :       gel(M, n+1) = gsub(gel(M, n+1),
     533       15484 :                          gmul(gel(M, n), gsub(A, uutoQ(n-1, d))));
     534         252 :   stripzeros(M); return M;
     535             : }
     536             : 
     537             : INLINE int
     538       14731 : RgV_is_CV(GEN x)
     539             : {
     540             :   long i;
     541       55147 :   for (i = lg(x)-1; i > 0; i--)
     542             :   {
     543       54559 :     long t = typ(gel(x,i));
     544       54559 :     if (!is_real_t(t) && t!= t_COMPLEX) return 0;
     545             :   }
     546         588 :   return 1;
     547             : }
     548             : 
     549             : static GEN
     550       14759 : get_Vga(GEN x, GEN *ldata)
     551             : {
     552       14759 :   if (typ(x)==t_VEC && RgV_is_CV(x)) { *ldata = NULL; return x; }
     553       14171 :   *ldata = lfunmisc_to_ldata_shallow_i(x);
     554       14171 :   if (*ldata) x = ldata_get_gammavec(*ldata);
     555       14171 :   return x;
     556             : }
     557             : GEN
     558          28 : gammamellininvasymp(GEN Vga, long nlim, long m)
     559             : {
     560          28 :   pari_sp av = avma;
     561             :   long status;
     562             :   GEN ldata;
     563          28 :   Vga = get_Vga(Vga, &ldata);
     564          28 :   if (!is_vec_t(typ(Vga)) || lg(Vga) == 1)
     565           7 :     pari_err_TYPE("gammamellininvasymp",Vga);
     566          21 :   return gerepilecopy(av, gammamellininvasymp_i(Vga, nlim, m, &status, 0));
     567             : }
     568             : 
     569             : /* Does the continued fraction of the asymptotic expansion M at oo of inverse
     570             :  * Mellin transform attached to Vga have zero Hankel determinants ? */
     571             : static long
     572        2891 : ishankelspec(GEN Vga)
     573             : {
     574        2891 :   long status, i, d = lg(Vga)-1;
     575             :   GEN M;
     576             : 
     577        2891 :   if (d == 5 || d == 7)
     578          21 :   { /* known bad cases: a x 5 and a x 7 */
     579         133 :     GEN v1 = gel(Vga, 1);
     580         588 :     for (i = 2; i <= d; ++i)
     581         476 :       if (!gequal(gel(Vga,i), v1)) break;
     582         133 :     if (i > d) return 1;
     583             :   }
     584        2758 :   else if (d==10 || d==14)
     585             :   { /* [ a x 5 , (a+1) x 5 ] */
     586           7 :     long d2 = d>>1;
     587           7 :     long s0 = 1, s1 = 0, sm1 = 0;
     588           7 :     GEN v1 = gel(Vga, 1);
     589          70 :     for (i = 2; i <= d; i++)
     590             :     {
     591          63 :       GEN s = gsub(gel(Vga,i),v1);
     592          63 :       if (gequal0(s)) s0++;
     593          35 :       else if(gequal1(s)) s1++;
     594           0 :       else if(gequalm1(s)) sm1++;
     595             :     }
     596           7 :     if (s0==d2 && (s1==d2 || sm1==d2)) return 1;
     597             :   }
     598             :   /* Heuristic: if 6 first terms in contfracinit don't fail, assume OK */
     599        2772 :   M = Klargeinit(Vga, 7, &status, 0);
     600        2772 :   return !contfracinit_i(M, 6);
     601             : }
     602             : 
     603             : /* Initialize data for computing m-th derivative of inverse Mellin */
     604             : GEN
     605       14738 : gammamellininvinit(GEN Vga, long m, long bitprec)
     606             : {
     607       14738 :   const double C2 = MELLININV_CUTOFF;
     608       14738 :   pari_sp ltop = avma;
     609             :   GEN A2, M, VS, VL, cd, ldata;
     610       14738 :   long nlimmax, status, d, prec = nbits2prec((4*bitprec)/3);
     611       14738 :   double E = M_LN2*bitprec, tmax = get_tmax(bitprec); /* = E/C2 */
     612             : 
     613       14738 :   if (m < 0)
     614           7 :     pari_err_DOMAIN("gammamellininvinit", "derivation order", "<", gen_0, stoi(m));
     615       14731 :   Vga = get_Vga(Vga, &ldata); d = lg(Vga)-1;
     616       14731 :   if (!is_vec_t(typ(Vga)) || !d) pari_err_TYPE("gammamellininvinit",Vga);
     617       14724 :   nlimmax = ceil(E * log2(1+M_PI*tmax) * C2 / get_D(d));
     618       14724 :   A2 = gaddsg(m*(2-d) + 1-d, sumVga(Vga));
     619       14724 :   cd = (d <= 2)? gen_2: gsqrt(gdivgu(int2n(d+1), d), nbits2prec(bitprec));
     620             :   /* if in Klarge, we have |t| > tmax = E/C2, thus nlim < E*C2/D. */
     621       14724 :   M = gammamellininvasymp_i(Vga, nlimmax, m, &status, prec);
     622       14724 :   if (status == 2)
     623             :   {
     624       11826 :     tmax = -1.; /* only use Klarge */
     625       11826 :     VS = gen_0;
     626             :   }
     627             :   else
     628             :   {
     629        2898 :     VS = Kderivsmallinit(ldata, Vga, m, bitprec);
     630        2898 :     if (status == 0 && ishankelspec(Vga)) status = 1;
     631        2898 :     if (status == 1)
     632             :     { /* a Hankel determinant vanishes => contfracinit is undefined.
     633             :          So compute K(t) / (1 - 1/(pi^2*t)) instead of K(t)*/
     634         133 :       GEN t = ginv(mppi(prec));
     635             :       long i;
     636       20902 :       for (i = 2; i < lg(M); ++i)
     637       20769 :         gel(M, i) = gadd(gel(M, i), gmul(gel(M, i-1), t));
     638             :     }
     639             :     else
     640        2765 :       M = RgC_gtofp(M, prec); /* convert from rationals to t_REAL: faster */
     641        2898 :     M = contfracinit(M, lg(M)-2);
     642             :   }
     643       14724 :   VL = mkvec3(mkvec2(M, stoi(status)), cd, A2);
     644       14724 :   return gerepilecopy(ltop, mkvec5(dbltor(tmax), Vga, stoi(m), VS, VL));
     645             : }
     646             : 
     647             : /* Compute m-th derivative of inverse Mellin at s2d = s^(d/2) using
     648             :  * initialization data. Use Taylor expansion at 0 for |s2d| < tmax, and
     649             :  * asymptotic expansion at oo otherwise. WARNING: assume that accuracy
     650             :  * has been increased according to tmax by the CALLING program. */
     651             : static GEN
     652     5638357 : gammamellininvrt_i(GEN K, GEN s, GEN s2d, long bit)
     653             : {
     654     5638357 :   if (dblmodulus(s2d) < GMi_get_tmax(K, bit))
     655       73971 :     return Kderivsmall(K, s, s2d, bit);
     656             :   else
     657     5564141 :     return Kderivlarge(K, s, s2d, bit);
     658             : }
     659             : GEN
     660     5632528 : gammamellininvrt(GEN K, GEN s2d, long bit)
     661     5632528 : { return gammamellininvrt_i(K, NULL, s2d, bit); }
     662             : 
     663             : /* Compute inverse Mellin at s. K from gammamellininv OR a Vga, in which
     664             :  * case the initialization data is computed. */
     665             : GEN
     666        5824 : gammamellininv(GEN K, GEN s, long m, long bitprec)
     667             : {
     668        5824 :   pari_sp av = avma;
     669             :   GEN s2d;
     670             :   long d;
     671             : 
     672        5824 :   if (!is_vec_t(typ(K)) || lg(K) != 6 || !is_vec_t(typ(GMi_get_Vga(K))))
     673          98 :     K = gammamellininvinit(K, m, bitprec);
     674        5824 :   d = GMi_get_degree(K);
     675        5824 :   s2d = gpow(s, gdivgu(gen_2, d), nbits2prec(bitprec));
     676        5824 :   return gerepileupto(av, gammamellininvrt_i(K, s, s2d, bitprec));
     677             : }

Generated by: LCOV version 1.14