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 23222-06b1652be) Lines: 311 323 96.3 %
Date: 2018-11-15 05:40:53 Functions: 30 31 96.8 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2015  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : 
      17             : /*******************************************************************/
      18             : /*               Computation of inverse Mellin                     */
      19             : /*               transforms of gamma products.                     */
      20             : /*******************************************************************/
      21             : #ifndef M_E
      22             : #define M_E 2.7182818284590452354
      23             : #endif
      24             : 
      25             : /* Handle complex Vga whose sum is real */
      26             : static GEN
      27       71935 : sumVga(GEN Vga) { return real_i(vecsum(Vga)); }
      28             : 
      29             : /* rough approximation to W0(a > -1/e), < 1% relative error */
      30             : double
      31       70755 : dbllambertW0(double a)
      32             : {
      33       70755 :   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       70755 :     double Wd = log(1.+a);
      43       70755 :     Wd *= (1.-log(Wd/a))/(1.+Wd);
      44       70755 :     if (a < 0.6482 && a > -0.1838) return Wd;
      45       70755 :     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      501899 : lemma526_i(double ac, double c, double t, double B)
      73             : {
      74      501899 :   double D = -B/ac; /* sgn(t) = sgn(a) = - sgn(D) */
      75      501899 :   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       53549 :     if (D < 100)
      93        9026 :       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       53549 :     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     1517991 : dblcoro526(double a, double c, double B)
     114             : {
     115     1517991 :   if (!a) return B <= 0? 0: pow(B/(2*M_PI*c), c);
     116      501899 :   if (B < 0) B = 1e-9;
     117      501899 :   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        5586 : MOD2(GEN x) { GEN q = gdivent(real_i(x),gen_2); return gsub(x,gmul2n(q,1)); }
     124             : static GEN
     125        1834 : RgV_MOD2(GEN v)
     126             : {
     127             :   long i, l;
     128        1834 :   GEN w = cgetg_copy(v,&l);
     129        1834 :   for (i=1; i<l; i++) gel(w,i) = MOD2(gel(v,i));
     130        1834 :   return w;
     131             : }
     132             : 
     133             : /* poles of the gamma factor */
     134             : static GEN
     135        1834 : gammapoles(GEN Vga)
     136             : {
     137        1834 :   long i, m, l = lg(Vga);
     138        1834 :   GEN P, B = RgV_MOD2(Vga), V = cgetg(l, t_VEC);
     139        1834 :   P = gen_indexsort(B, (void*)lexcmp, cmp_nodata);
     140        6517 :   for (i = m = 1; i < l;)
     141             :   {
     142        2849 :     GEN u = gel(B, P[i]);
     143             :     long k;
     144        5586 :     for(k = i+1; k < l; k++)
     145        3752 :       if (gexpo(gsub(u, gel(B, P[k]))) > -32) break;
     146        2849 :     gel(V, m++) = vecpermute(Vga, vecslice(P,i,k-1));
     147        2849 :     i = k;
     148             :   }
     149        1834 :   setlg(V, m); return V;
     150             : }
     151             : 
     152             : static GEN
     153      186828 : sercoeff(GEN x, long n, long prec)
     154             : {
     155      186828 :   long N = n - valp(x);
     156      186828 :   return (N < 0)? gen_0: gtofp(gel(x, N+2), prec);
     157             : }
     158             : 
     159             : /* generalized power series expansion of inverse Mellin around x = 0;
     160             :  * m-th derivative */
     161             : static GEN
     162        1834 : Kderivsmallinit(GEN Vga, long m, long bitprec)
     163             : {
     164        1834 :   double C2 = MELLININV_CUTOFF;
     165        1834 :   long N, j, l, d = lg(Vga)-1, limn, prec = nbits2prec(1+bitprec*(1+M_PI*d/C2));
     166             :   GEN LA, L, M, mat;
     167             : 
     168        1834 :   Vga = gprec_wensure(Vga, prec);
     169        1834 :   LA = gammapoles(Vga); N = lg(LA)-1;
     170        1834 :   L = cgetg(N+1, t_VECSMALL);
     171        1834 :   M = cgetg(N+1, t_VEC);
     172        4683 :   for (j = 1; j <= N; ++j)
     173             :   {
     174        2849 :     GEN S = gel(LA,j);
     175        2849 :     L[j] = lg(S)-1;
     176        2849 :     gel(M,j) = gsubsg(2, gel(S, vecindexmin(real_i(S))));
     177             :   }
     178        1834 :   limn = ceil(2*M_LN2*bitprec/(d*dbllambertW0(C2/(M_PI*M_E))));
     179        1834 :   mat = cgetg(N+1, t_VEC);
     180        1834 :   l = limn + 2;
     181        4683 :   for (j=1; j <= N; j++)
     182             :   {
     183        2849 :     GEN C, c, mj = gel(M,j), pr = gen_1, t = gen_1;
     184        2849 :     long i, k, n, lj = L[j], lprecdl = lj+3;
     185       12117 :     for (i = 1; i <= d; i++)
     186             :     {
     187        9268 :       GEN a = gmul2n(gadd(mj, gel(Vga,i)), -1);
     188        9268 :       GEN u = deg1pol_shallow(ghalf, a, 0);
     189        9268 :       pr = gmul(pr, ggamma(RgX_to_ser(u, lprecdl), prec));
     190        9268 :       t = gmul(t, u);
     191             :     }
     192        2849 :     c = cgetg(limn+2,t_COL); gel(c,1) = pr;
     193      101816 :     for (n=1; n <= limn; n++)
     194             :     {
     195       98967 :       GEN T = RgX_translate(t, stoi(-2*n));
     196       98967 :       if (n == 1) gel(T,2) = gen_0; /* in case t inexact */
     197       98967 :       gel(c,n+1) = gdiv(gel(c,n), T);
     198             :     }
     199             : 
     200        2849 :     gel(mat, j) = C = cgetg(lj+1, t_COL);
     201        8435 :     for (k = 1; k <= lj; k++)
     202             :     {
     203        5586 :       GEN L = cgetg(l, t_POL);
     204        5586 :       for (n = 2; n < l; n++) gel(L,n) = sercoeff(gel(c,n), -k, prec);
     205        5586 :       L[1] = evalsigne(1)|evalvarn(0); gel(C,k) = L;
     206             :     }
     207             :     /* C[k] = \sum_n c_{j,k} t^n =: C_k(t) in Dokchitser's Algo 3.3 */
     208             :     /* Take m-th derivative of t^(-M+2) sum_k (-ln t)^k/k! C_k(t^2) */
     209        2849 :     if (m)
     210             :     {
     211         119 :       mj = gsubgs(mj, 2);
     212         238 :       for (i = 1; i <= m; i++, mj = gaddgs(mj,1))
     213         322 :         for (k = 1; k <= lj; k++)
     214             :         {
     215         203 :           GEN c = gel(C,k), d = RgX_shift_shallow(gmul2n(RgX_deriv(c),1), 1);
     216         203 :           c = RgX_Rg_mul(c, mj);
     217         203 :           if (k < lj) c = RgX_add(c, gel(C,k+1));
     218         203 :           gel(C,k) = RgX_sub(d, c);
     219             :         }
     220         119 :       gel(M,j) = gaddgs(mj,2);
     221             :     }
     222        8435 :     for (k = 1; k <= lj; k++)
     223             :     {
     224        5586 :       GEN c = gel(C,k);
     225        5586 :       if (k > 2) c = RgX_Rg_div(c, mpfact(k-1));
     226        5586 :       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        1834 :   return mkvec3(L, RgV_neg(M), 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      208062 : evalvec(GEN vec, long lim, GEN u, GEN ui)
     237             : {
     238      208062 :   pari_sp ltop = avma;
     239      208062 :   GEN S = gen_0;
     240             :   long n;
     241      208062 :   lim = minss(lim, lg(vec)-1);
     242      208062 :   if (!ui)
     243       54105 :     for (n = lim; n >= 1; --n) S = gmul(u, gadd(gel(vec,n), S));
     244             :   else
     245             :   {
     246      153957 :     for (n = 1; n <= lim; ++n) S = gmul(ui, gadd(gel(vec,n), S));
     247      153957 :     S = gmul(gpowgs(u, n), S);
     248             :   }
     249      208062 :   return gerepileupto(ltop, S);
     250             : }
     251             : 
     252             : /* gammamellininvinit accessors */
     253             : static double
     254     4598553 : get_tmax(long bitprec)
     255     4598553 : { return (M_LN2 / MELLININV_CUTOFF) * bitprec ; }
     256             : static GEN
     257     9917977 : GMi_get_Vga(GEN K) { return gel(K,2); }
     258             : static long
     259     4983210 : GMi_get_m(GEN K) { return itos( gel(K,3) ); }
     260             : static GEN /* [lj,mj,mat], Kderivsmall only */
     261     5043105 : GMi_get_VS(GEN K) { return gel(K,4); }
     262             : static GEN /* [Ms,cd,A2], Kderivlarge only */
     263     9846630 : GMi_get_VL(GEN K) { return gel(K,5); }
     264             : static double
     265     4983210 : GMi_get_tmax(GEN K, long bitprec)
     266     4983210 : { 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 = M_LN2*bitprec/MELLININV_CUTOFF*/
     271             : static GEN
     272       59895 : Kderivsmall(GEN K, GEN x, GEN x2d, long bitprec)
     273             : {
     274       59895 :   pari_sp ltop = avma;
     275       59895 :   GEN Vga = GMi_get_Vga(K), VS = GMi_get_VS(K);
     276       59895 :   GEN L = gel(VS,1), M = gel(VS,2), mat = gel(VS,3);
     277             :   GEN d2, Lx, x2, x2i, A, S, pi;
     278       59895 :   long prec, j, k, limn, m = GMi_get_m(K), N = lg(L)-1, d = lg(Vga)-1;
     279       59895 :   double xd, Wd, Ed = M_LN2*bitprec / d;
     280             : 
     281       59895 :   xd = maxdd(M_PI*dblmodulus(x2d), 1E-13); /* pi |x|^2/d unless x tiny */
     282       59895 :   if (xd > Ed) pari_err_BUG("Kderivsmall (x2d too large)");
     283             :   /* Lemma 5.2.6 (2), a = 1 + log(Pi x^(2/d)) = log(e / xd),
     284             :    * B = log(2)*bitprec / d = Ed */
     285       59895 :   Wd = dbllambertW0( Ed / (M_E*xd) ); /* solution of w exp(w) = B exp(-a)*/
     286       59895 :   limn = (long) ceil(2*Ed/Wd);
     287       59895 :   prec = nbits2prec((long) ceil(bitprec+d*xd/M_LN2));
     288       59895 :   pi = mppi(prec);
     289       59895 :   d2 = gdivsg(d,gen_2);
     290       59895 :   if (x)
     291        1052 :     x = gmul(gtofp(x,prec), gpow(pi,d2,prec));
     292             :   else
     293       58843 :     x = gpow(gmul(gtofp(x2d,prec),pi), d2, prec);
     294             :   /* at this stage, x has been replaced by pi^(d/2) x */
     295       59895 :   x2 = gsqr(x);
     296       59895 :   Lx = gpowers(gneg(glog(x,prec)), vecsmall_max(L));
     297       59895 :   x2i = (gcmp(gnorml2(x2), gen_1) <= 0)? NULL: ginv(x2);
     298       59895 :   S = gen_0;
     299      156195 :   for (j = 1; j <= N; ++j)
     300             :   {
     301       96300 :     long lj = L[j];
     302       96300 :     GEN s = gen_0;
     303      304362 :     for (k = 1; k <= lj; k++)
     304      208062 :       s = gadd(s, gmul(gel(Lx,k), evalvec(gmael(mat,j,k), limn, x2, x2i)));
     305       96300 :     S = gadd(S, gmul(gpow(x, gel(M,j), prec), s));
     306             :   }
     307       59895 :   A = gsubsg(m*d, sumVga(Vga));
     308       59895 :   if (!gequal0(A)) S = gmul(S, gpow(pi, gmul2n(A,-1), prec));
     309       59895 :   return gerepileupto(ltop, gtofp(S, nbits2prec(bitprec)));
     310             : }
     311             : 
     312             : /* In Klarge, we conpute K(t) as (asymptotic) * F(z), where F ~ 1 is given by
     313             :  * a continued fraction and z = Pi t^(2/d). If we take 2n terms in F (n terms
     314             :  * in Euler form), F_n(z) - F(z) is experimentally in exp(- C sqrt(n*z))
     315             :  * where C ~ 8 for d > 2 [HEURISTIC] and C = 4 (theorem) for d = 1 or d = 2
     316             :  * and vga = [0,1]. For e^(-E) absolute error, we want
     317             :  *   exp(-C sqrt(nz)) < e^-(E+a), where a ~ ln(asymptotic)
     318             :  * i.e.  2n > (E+a)^2 / t^(2/d) * 2/(C^2 Pi); C^2*Pi/2 ~ 100.5 ~ 101
     319             :  *
     320             :  * In fact, this model becomes wrong for z large: we use instead
     321             :  *
     322             :  *   exp(- sqrt(D * nz/log(z+1))) < e^-(E+a),
     323             :  * i.e.  2n > (E+a)^2 * log(1 + Pi t^(2/d))/ t^(2/d) * 2/(D Pi); */
     324             : static double
     325     4935103 : get_D(long d) { return d <= 2 ? 157. : 180.; }
     326             : /* if (abs), absolute error rather than relative */
     327             : static void
     328     4923315 : Kderivlarge_optim(GEN K, long abs, GEN t2d,GEN gcd, long *pbitprec, long *pnlim)
     329             : {
     330     4923315 :   GEN Vga = GMi_get_Vga(K), VL = GMi_get_VL(K), A2 = gel(VL,3);
     331     4923315 :   long bitprec = *pbitprec, d = lg(Vga)-1;
     332     4923315 :   const double D = get_D(d), td = dblmodulus(t2d), cd = gtodouble(gcd);
     333     4923315 :   double a, rtd, E = M_LN2*bitprec;
     334             : 
     335     4923315 :   rtd = (typ(t2d) == t_COMPLEX)? gtodouble(gel(t2d,1)): td;
     336             : 
     337             :   /* A2/2 = A, log(td) = (2/d)*log t */
     338     4923315 :   a = d*gtodouble(A2)*log2(td)/2 - (M_PI/M_LN2)*d*rtd + log2(cd); /*log2 K(t)~a*/
     339             : 
     340             :   /* if bitprec <= 0, caller should return K(t) ~ 0 */
     341     4923315 :   bitprec += 64;
     342     4923315 :   if (abs)
     343             :   {
     344     4918634 :     bitprec += ceil(a);
     345     4918634 :     if (a <= -65) E = M_LN2*bitprec; /* guarantees E <= initial E */
     346             :   }
     347     4923315 :   *pbitprec = bitprec;
     348     4923315 :   *pnlim = ceil(E*E * log2(1+M_PI*td) / (D*td));
     349     4923315 : }
     350             : 
     351             : /* Compute m-th derivative of inverse Mellin at t by continued fraction of
     352             :  * asymptotic expansion; t2d = t^(2/d). If t is NULL, "lfun" mode: don't
     353             :  * bother about complex branches + use absolute (rather than relative)
     354             :  * accuracy */
     355             : static GEN
     356     4923315 : Kderivlarge(GEN K, GEN t, GEN t2d, long bitprec0)
     357             : {
     358     4923315 :   pari_sp ltop = avma;
     359     4923315 :   GEN tdA, P, S, pi, z, Vga = GMi_get_Vga(K);
     360     4923315 :   const long d = lg(Vga)-1;
     361     4923315 :   GEN M, VL = GMi_get_VL(K), Ms = gel(VL,1), cd = gel(VL,2), A2 = gel(VL,3);
     362     4923315 :   long status, prec, nlim, m = GMi_get_m(K), bitprec = bitprec0;
     363             : 
     364     4923315 :   Kderivlarge_optim(K, !t, t2d, cd, &bitprec, &nlim);
     365     4923315 :   if (bitprec <= 0) return gen_0;
     366     4856122 :   prec = nbits2prec(bitprec);
     367     4856122 :   t2d = gtofp(t2d, prec);
     368     4856122 :   if (t)
     369        4681 :     tdA = gpow(t, gdivgs(A2,d), prec);
     370             :   else
     371     4851441 :     tdA = gpow(t2d, gdivgs(A2,2), prec);
     372     4856122 :   tdA = gmul(cd, tdA);
     373             : 
     374     4856122 :   pi = mppi(prec);
     375     4856122 :   z = gmul(pi, t2d);
     376     4856122 :   P = gmul(tdA, gexp(gmulsg(-d, z), prec));
     377     4856122 :   if (m) P = gmul(P, gpowgs(mulsr(-2, pi), m));
     378     4856122 :   M = gel(Ms,1);
     379     4856122 :   status = itos(gel(Ms,2));
     380     4856122 :   if (status == 2)
     381             :   {
     382      394303 :     if (lg(M) == 2) /* shortcut: continued fraction is constant */
     383      393974 :       S = gel(M,1);
     384             :     else
     385         329 :       S = poleval(RgV_to_RgX(M, 0), ginv(z));
     386             :   }
     387             :   else
     388             :   {
     389     4461819 :     S = contfraceval_inv(M, z, nlim/2);
     390     4461819 :     if (DEBUGLEVEL>3)
     391             :     {
     392           0 :       GEN S0 = contfraceval_inv(M, z, nlim/2 + 1);
     393           0 :       long e = gexpo(gmul(P, gabs(gsub(S,S0),0)));
     394           0 :       if (-e < bitprec0)
     395           0 :         err_printf("Kderivlarge: e = %ld, bit = %ld\n",e,bitprec0);
     396             :     }
     397     4461819 :     if (status == 1) S = gmul(S, gsubsg(1, ginv(gmul(z, pi))));
     398             :   }
     399     4856122 :   return gerepileupto(ltop, gmul(P, S));
     400             : }
     401             : 
     402             : /* Dokchitser's coefficients used for asymptotic expansion of inverse Mellin
     403             :  * 2 <= p <= min(n+1, d) */
     404             : static GEN
     405      487371 : fun_vp(long p, long n, long d, GEN SM, GEN vsinh)
     406             : {
     407      487371 :   pari_sp ltop = avma;
     408             :   long m, j, k;
     409      487371 :   GEN s = gen_0;
     410     2306203 :   for (m = 0; m <= p; ++m)
     411             :   {
     412     1818832 :     GEN pr = gen_1, s2 = gen_0, sh = gel(vsinh, d-p+1);/* (sh(x)/x)^(d-p) */
     413     1818832 :     long pm = p-m;
     414     1818832 :     for (j = m; j < p; ++j) pr = muliu(pr, d-j);
     415     4606562 :     for (k = 0; k <= pm; k+=2)
     416             :     {
     417     2787730 :       GEN e = gdiv(powuu(2*n-p+1, pm-k), mpfact(pm-k));
     418     2787730 :       s2 = gadd(s2, gmul(e, RgX_coeff(sh, k)));
     419             :     }
     420     1818832 :     s = gadd(s, gmul(gmul(gel(SM, m+1), pr), s2));
     421     1818832 :     if (gc_needed(ltop, 1)) s = gerepilecopy(ltop, s);
     422             :   }
     423      487371 :   return gerepileupto(ltop, gmul(gdivsg(-d, powuu(2*d, p)), s));
     424             : }
     425             : 
     426             : /* Asymptotic expansion of inverse Mellin, to length nlimmax. Set status = 0
     427             :  * (regular), 1 (one Hankel determinant vanishes => contfracinit will fail)
     428             :  * or 2 (same as 1, but asymptotic expansion is finite!)
     429             :  *
     430             :  * If status = 2, the asymptotic expansion is finite so return only
     431             :  * the necessary number of terms nlim <= nlimmax + d. */
     432             : static GEN
     433       11802 : Klargeinit0(GEN Vga, long nlimmax, long *status)
     434             : {
     435       11802 :   const long prec = LOWDEFAULTPREC;
     436       11802 :   const long d = lg(Vga)-1;
     437             :   long k, n, m, cnt;
     438             :   GEN pol, SM, nS1, se, vsinh, M, dk;
     439             : 
     440       11802 :   if (d==1 || (d==2 && gequal1(gabs(gsub(gel(Vga,1), gel(Vga,2)), prec))))
     441             :   { /* shortcut */
     442        9933 :     *status = 2; return mkvec(gen_1);
     443             :   }
     444             :   /* d >= 2 */
     445        1869 :   *status = 0;
     446        1869 :   pol = roots_to_pol(gneg(Vga), 0); /* deg(pol) = d */
     447        1869 :   nS1 = gpowers(gneg(RgX_coeff(pol, d-1)), d);
     448        1869 :   dk = gpowers(utoi(d), d-1);
     449        1869 :   SM = cgetg(d+3, t_VEC);
     450        9408 :   for (m = 0; m <= d; ++m)
     451             :   {
     452        7539 :     pari_sp btop = avma;
     453        7539 :     GEN s = gmul(gdivgs(gel(nS1, m+1), d), binomialuu(d, m));
     454       20048 :     for (k = 1; k <= m; ++k)
     455             :     {
     456       12509 :       GEN e = gmul(gel(nS1, m-k+1), gel(dk, k));
     457       12509 :       s = gadd(s, gmul(gmul(e, binomialuu(d-k, m-k)), RgX_coeff(pol, d-k)));
     458             :     }
     459        7539 :     gel(SM, m+1) = gerepileupto(btop, s);
     460             :   }
     461        1869 :   se = gdiv(gsinh(RgX_to_ser(pol_x(0), d+2), prec), pol_x(0));
     462        1869 :   vsinh = gpowers(se, d);
     463        1869 :   M = vectrunc_init(nlimmax + d);
     464        1869 :   vectrunc_append(M, gen_1);
     465      248558 :   for (n=2, cnt=0; (n <= nlimmax) || cnt; ++n)
     466             :   {
     467      248558 :     pari_sp btop = avma;
     468      248558 :     long p, ld = minss(d, n);
     469      248558 :     GEN s = gen_0;
     470      735929 :     for (p = 2; p <= ld; ++p)
     471      487371 :       s = gadd(s, gmul(fun_vp(p, n-1, d, SM, vsinh), gel(M, n+1-p)));
     472      248558 :     s = gerepileupto(btop, gdivgs(s, n-1));
     473      248558 :     vectrunc_append(M, s);
     474      248558 :     if (!isintzero(s))
     475             :     {
     476      248509 :       if (n >= nlimmax) break;
     477      246668 :       cnt = 0;
     478             :     }
     479             :     else
     480             :     {
     481          49 :       cnt++; *status = 1;
     482          49 :       if (cnt >= d-1) { *status = 2; setlg(M, lg(M) - (d-1)); break; }
     483             :     }
     484             :   }
     485        1869 :   return M;
     486             : }
     487             : 
     488             : /* remove trailing zeros from vector. */
     489             : static void
     490         252 : stripzeros(GEN M)
     491             : {
     492             :   long i;
     493         441 :   for(i = lg(M)-1; i >= 1; --i)
     494         441 :     if (!gequal0(gel(M, i))) break;
     495         252 :   setlg(M, i+1);
     496         252 : }
     497             : 
     498             : /* Asymptotic expansion of the m-th derivative of inverse Mellin, to length
     499             :  * nlimmax. If status = 2, the asymptotic expansion is finite so return only
     500             :  * the necessary number of terms nlim <= nlimmax + d. */
     501             : static GEN
     502       11802 : gammamellininvasymp_i(GEN Vga, long nlimmax, long m, long *status)
     503             : {
     504       11802 :   pari_sp ltop = avma;
     505             :   GEN M, A, Aadd;
     506             :   long d, i, nlim, n;
     507             : 
     508       11802 :   M = Klargeinit0(Vga, nlimmax, status);
     509       11802 :   if (!m) return gerepilecopy(ltop, M);
     510         252 :   d = lg(Vga)-1;
     511             :   /* half the exponent of t in asymptotic expansion. */
     512         252 :   A = gdivgs(gaddsg(1-d, sumVga(Vga)), 2*d);
     513         252 :   if (*status == 2) M = shallowconcat(M, zerovec(m));
     514         252 :   nlim = lg(M)-1;
     515         252 :   Aadd = gdivgs(stoi(2-d), 2*d); /* (1/d) - (1/2) */
     516         532 :   for (i = 1; i <= m; i++, A = gadd(A,Aadd))
     517        8022 :     for (n = nlim-1; n >= 1; --n)
     518       15484 :       gel(M, n+1) = gsub(gel(M, n+1),
     519        7742 :                          gmul(gel(M, n), gsub(A, gdivgs(stoi(n-1), d))));
     520         252 :   stripzeros(M);
     521         252 :   return gerepilecopy(ltop, M);
     522             : }
     523             : GEN
     524          14 : gammamellininvasymp(GEN Vga, long nlimmax, long m)
     525             : { long status;
     526          14 :   if (!is_vec_t(typ(Vga))) pari_err_TYPE("gammamellininvinit",Vga);
     527          14 :   return gammamellininvasymp_i(Vga, nlimmax, m, &status);
     528             : }
     529             : 
     530             : /* Does the continued fraction of the asymptotic expansion M at oo of inverse
     531             :  * Mellin transform attached to Vga have zero Hankel determinants ? */
     532             : static long
     533        1827 : ishankelspec(GEN Vga, GEN M)
     534             : {
     535        1827 :   long status, i, d = lg(Vga)-1;
     536             : 
     537        1827 :   if (d == 5 || d == 7)
     538             :   { /* known bad cases: a x 5 and a x 7 */
     539          70 :     GEN v1 = gel(Vga, 1);
     540         294 :     for (i = 2; i <= d; ++i)
     541         238 :       if (!gequal(gel(Vga,i), v1)) break;
     542          70 :     if (i > d) return 1;
     543             :   }
     544        1771 :   status = 0;
     545             :   /* Heuristic: if 6 first terms in contfracinit don't fail, assume it's OK */
     546        1771 :   pari_CATCH(e_INV) { status = 1; }
     547        1771 :   pari_TRY { contfracinit(M, minss(lg(M)-2,6)); }
     548        1771 :   pari_ENDCATCH; return status;
     549             : }
     550             : 
     551             : /* Initialize data for computing m-th derivative of inverse Mellin */
     552             : GEN
     553       11788 : gammamellininvinit(GEN Vga, long m, long bitprec)
     554             : {
     555       11788 :   pari_sp ltop = avma;
     556             :   GEN A2, M, VS, VL, cd;
     557       11788 :   long d = lg(Vga)-1, status;
     558       11788 :   const double C2 = MELLININV_CUTOFF, D = get_D(d);
     559       11788 :   double E = M_LN2*bitprec, tmax = get_tmax(bitprec); /* = E/C2 */
     560       11788 :   const long nlimmax = ceil(E*log2(1+M_PI*tmax)*C2/D);
     561             : 
     562       11788 :   if (!is_vec_t(typ(Vga))) pari_err_TYPE("gammamellininvinit",Vga);
     563       11788 :   A2 = gaddsg(m*(2-d) + 1-d, sumVga(Vga));
     564       11788 :   cd = (d <= 2)? gen_2: gsqrt(gdivgs(int2n(d+1), d), nbits2prec(bitprec));
     565             :   /* if in Klarge, we have |t| > tmax = E/C2, thus nlim < E*C2/D. */
     566       11788 :   M = gammamellininvasymp_i(Vga, nlimmax, m, &status);
     567       11788 :   if (status == 2)
     568             :   {
     569        9954 :     tmax = -1.; /* only use Klarge */
     570        9954 :     VS = gen_0;
     571             :   }
     572             :   else
     573             :   {
     574        1834 :     long prec = nbits2prec((4*bitprec)/3);
     575        1834 :     VS = Kderivsmallinit(Vga, m, bitprec);
     576        1834 :     if (status == 0 && ishankelspec(Vga, M)) status = 1;
     577        1834 :     if (status == 1)
     578             :     { /* a Hankel determinant vanishes => contfracinit is undefined.
     579             :          So compute K(t) / (1 - 1/(pi^2*t)) instead of K(t)*/
     580          63 :       GEN t = ginv(mppi(prec));
     581             :       long i;
     582        9713 :       for (i = 2; i < lg(M); ++i)
     583        9650 :         gel(M, i) = gadd(gel(M, i), gmul(gel(M, i-1), t));
     584             :     }
     585             :     else
     586        1771 :       M = RgC_gtofp(M, prec); /* convert from rationals to t_REAL: faster */
     587        1834 :     M = contfracinit(M, lg(M)-2);
     588             :   }
     589       11788 :   VL = mkvec3(mkvec2(M, stoi(status)), cd, A2);
     590       11788 :   return gerepilecopy(ltop, mkvec5(dbltor(tmax), Vga, stoi(m), VS, VL));
     591             : }
     592             : 
     593             : /* Compute m-th derivative of inverse Mellin at s2d = s^(d/2) using
     594             :  * initialization data. Use Taylor expansion at 0 for |s2d| < tmax, and
     595             :  * asymptotic expansion at oo otherwise. WARNING: assume that accuracy
     596             :  * has been increased according to tmax by the CALLING program. */
     597             : GEN
     598     4977477 : gammamellininvrt(GEN K, GEN s2d, long bitprec)
     599             : {
     600     4977477 :   if (dblmodulus(s2d) < GMi_get_tmax(K, bitprec))
     601       58843 :     return Kderivsmall(K, NULL, s2d, bitprec);
     602             :   else
     603     4918634 :     return Kderivlarge(K, NULL, s2d, bitprec);
     604             : }
     605             : 
     606             : /* Compute inverse Mellin at s. K from gammamellininv OR a Vga, in which
     607             :  * case the initialization data is computed. */
     608             : GEN
     609        5733 : gammamellininv(GEN K, GEN s, long m, long bitprec)
     610             : {
     611        5733 :   pari_sp av = avma;
     612             :   GEN z, s2d;
     613             :   long d;
     614        5733 :   if (!is_vec_t(typ(K))) pari_err_TYPE("gammamellininv",K);
     615        5733 :   if (lg(K) != 6 || !is_vec_t(typ(GMi_get_Vga(K))))
     616          14 :     K = gammamellininvinit(K, m, bitprec);
     617        5733 :   d = lg(GMi_get_Vga(K))-1;
     618        5733 :   s2d = gpow(s, gdivgs(gen_2, d), nbits2prec(bitprec));
     619        5733 :   if (dblmodulus(s2d) < GMi_get_tmax(K, bitprec))
     620        1052 :     z = Kderivsmall(K, s, s2d, bitprec);
     621             :   else
     622        4681 :     z = Kderivlarge(K, s, s2d, bitprec);
     623        5733 :   return gerepileupto(av, z);
     624             : }

Generated by: LCOV version 1.13