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 - hypergeom.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.14.0 lcov report (development 27775-aca467eab2) Lines: 666 675 98.7 %
Date: 2022-07-03 07:33:15 Functions: 70 70 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2017  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             : /********************************************************************/
      16             : /**                                                                **/
      17             : /**                   HYPERGEOMETRIC FUNCTIONS                     **/
      18             : /**                                                                **/
      19             : /********************************************************************/
      20             : 
      21             : #include "pari.h"
      22             : #include "paripriv.h"
      23             : 
      24             : static GEN
      25         350 : F10(GEN a, GEN z, long prec)
      26         350 : { return gpow(gsubsg(1,z), gneg(a), prec); }
      27             : 
      28             : static int
      29      232974 : isnegint2(GEN a, long *pa)
      30             : {
      31             :   GEN b;
      32      232974 :   if (!gequal0(imag_i(a))) return 0;
      33      189462 :   a = real_i(a); if (gsigne(a) > 0) return 0;
      34       42560 :   b = ground(a); if (!gequal(a, b)) return 0;
      35       40712 :   if (pa) *pa = -itos(b);
      36       40712 :   return 1;
      37             : }
      38             : static int
      39      145096 : isnegint(GEN a) { return isnegint2(a, NULL); }
      40             : static int
      41       73696 : is0(GEN a, long bit) { return gequal0(a) || gexpo(a) < -bit; }
      42             : static int
      43       12978 : islong(GEN a, long *m, long prec)
      44             : {
      45       12978 :   *m = itos(ground(real_i(a)));
      46       12978 :   if (is0(gsubgs(a, *m), prec2nbits(prec))) return 1;
      47        4431 :   return 0;
      48             : }
      49             : 
      50             : static GEN
      51          21 : F01(GEN a, GEN z, long prec)
      52             : {
      53             :   GEN A, B, al, sz;
      54          21 :   if (is0(z, prec2nbits(prec)-5)) return real_1(prec);
      55          21 :   sz = gsqrt(z, prec); al = gsubgs(a, 1);
      56          21 :   A = gmul(ggamma(a, prec), gpow(sz, gneg(al), prec));
      57          21 :   B = ibessel(al, gmul2n(sz,1), prec);
      58          21 :   return isexactzero(imag_i(z))? mulreal(A,B): gmul(A,B);
      59             : }
      60             : 
      61             : /* Airy functions [Ai,Bi] */
      62             : static GEN
      63          70 : airy_i(GEN x, long prec)
      64             : {
      65          70 :   long bit = prec2nbits(prec), tx = typ(x), prec2;
      66             :   GEN a, b, A, B, z, z2;
      67          70 :   if (!is_scalar_t(tx)) pari_err_TYPE("airy",x);
      68          63 :   if (is0(x, bit))
      69             :   {
      70           7 :     GEN s = sqrtnr_abs(utor(3,prec), 6), s3 = powrs(s,3), s4 = mulrr(s,s3);
      71           7 :     A = invr(mulrr(s4, ggamma(uutoQ(2,3), prec)));
      72           7 :     B = mulrr(A, s3); return mkvec2(A, B);
      73             :   }
      74          56 :   prec2 = prec + EXTRAPRECWORD;
      75          56 :   x = gprec_wensure(x, prec2);
      76          56 :   z = gsqrt(gpowgs(x,3), prec2); z2 = gdivgu(gmul2n(z,1),3);
      77          56 :   if (is_real_t(tx) && gsigne(x) > 0)
      78          42 :     a = b = gsqrt(x, prec2); /* expression simplifies */
      79             :   else
      80             :   {
      81          14 :     a = gsqrtn(z, utoipos(3), NULL, prec2);
      82          14 :     b = gdiv(x, a);
      83             :   }
      84          56 :   a = gmul(a, ibessel(sstoQ(-1,3),z2, prec));
      85          56 :   b = gmul(b, ibessel(uutoQ(1,3), z2, prec));
      86          56 :   if (isexactzero(imag_i(x))) { a = real_i(a); b = real_i(b); }
      87          56 :   A = gdivgu(gsub(a,b), 3);
      88          56 :   B = gdiv(gadd(a,b), sqrtr_abs(utor(3, prec)));
      89             : 
      90          56 :   bit -= gexpo(a) + 16;
      91          56 :   if (!is0(A, bit) && !is0(B, bit)) return mkvec2(A, B);
      92          28 :   prec = precdbl(prec); x = gprec_wensure(x, prec); return airy_i(x, prec);
      93             : }
      94             : GEN
      95          42 : airy(GEN z, long prec)
      96          42 : { pari_sp av = avma; return gerepilecopy(av, airy_i(z, prec)); }
      97             : 
      98             : /* Gamma(a)*Gamma(b) */
      99             : static GEN
     100       28868 : mulgamma2(GEN a, GEN b, long prec)
     101       28868 : { return gmul(ggamma(a, prec), ggamma(b, prec)); }
     102             : /* Gamma(a)/Gamma(b) */
     103             : static GEN
     104         140 : divgamma2(GEN a, GEN b, long prec)
     105         140 : { return gdiv(ggamma(a, prec), ggamma(b, prec)); }
     106             : /* Gamma(v[1])*Gamma(v[2]) */
     107             : static GEN
     108       28546 : mulgammav2(GEN v, long prec)
     109       28546 : { return mulgamma2(gel(v,1), gel(v,2), prec); }
     110             : 
     111             : /***********************************************************************/
     112             : /**                 CONFLUENT HYPERGEOMETRIC U(a,b,z)                 **/
     113             : /***********************************************************************/
     114             : static GEN Ftaylor(GEN N, GEN D, GEN z, long prec);
     115             : /* b not integral; use 1F1; prec0 is precision we really want */
     116             : static GEN
     117          70 : hyperu_F11(GEN a, GEN b, GEN z, long prec0, long prec)
     118             : {
     119          70 :   GEN S1, S2, b1 = gsubsg(1, b), ab1 = gadd(a, b1);
     120          70 :   if (isnegint(ab1)) S1 = gen_0;
     121             :   else
     122             :   {
     123          70 :     GEN tmp = Ftaylor(mkvec(a), mkvec(b), z, prec);
     124          70 :     S1 = gmul(divgamma2(b1, ab1, prec), tmp);
     125             :   }
     126          70 :   if (isnegint(a)) S2 = gen_0;
     127             :   else
     128             :   {
     129          70 :     GEN tmp = Ftaylor(mkvec(ab1), mkvec(gaddsg(1, b1)), z, prec);
     130          70 :     S2 = gmul(divgamma2(gneg(b1), a, prec), tmp);
     131          70 :     S2 = gmul(S2, gpow(z, b1, prec));
     132             :   }
     133          70 :   S1 = gadd(S1, S2);
     134          70 :   if (gexpo(S1)-gexpo(S2) >= prec2nbits(prec0) - prec2nbits(prec))
     135          35 :     return S1;
     136          35 :   prec = precdbl(prec);
     137          35 :   a = gprec_wensure(a, prec);
     138          35 :   b = gprec_wensure(b, prec);
     139          35 :   z = gprec_wensure(z, prec);
     140          35 :   return hyperu_F11(a, b, z, prec0, prec);
     141             : }
     142             : /* one branch of this must assume x > 0 (a,b complex); see Temme, The
     143             :  * numerical computation of the confluent hypergeometric function U(a,b,z),
     144             :  * Numer. Math. 41 (1983), no. 1, 63-82. */
     145             : static GEN
     146          70 : hyperu_i(GEN a, GEN b, GEN x, long prec)
     147             : {
     148             :   GEN u, S, P, T, zf, a1;
     149             :   long k, n, bit, l, bigx;
     150             : 
     151          70 :   if (gequal0(imag_i(x))) x = real_i(x);
     152          70 :   l = precision(x); if (!l) l = prec;
     153          70 :   prec = l;
     154          70 :   a1 = gaddsg(1, gsub(a,b));
     155          70 :   P = gmul(a1, a);
     156          70 :   S = gadd(a1, a);
     157          70 :   n = (long)(prec2nbits_mul(l, M_LN2) + M_PI*sqrt(dblmodulus(P)));
     158          70 :   bigx = dbllog2(x) >= log2((double)n);
     159          70 :   if (!bigx && (!isint(b,&b) || typ(x) == t_COMPLEX || gsigne(x) <= 0))
     160             :   {
     161          35 :     if (typ(b) == t_INT)
     162             :     {
     163          14 :       bit = prec2nbits(l); l += l-2;
     164          14 :       b = gadd(b, real2n(-bit, l));
     165          14 :       a = gprec_wensure(a, l);
     166          14 :       x = gprec_wensure(x, l);
     167             :     }
     168          35 :     return hyperu_F11(a, b, x, l, l);
     169             :   }
     170          35 :   bit = prec2nbits(l)-1;
     171          35 :   l += EXTRAPRECWORD;
     172          35 :   T = gadd(gadd(P, gmulsg(n-1, S)), sqru(n-1));
     173          35 :   x = gtofp(x, l);
     174          35 :   if (!bigx)
     175             :   { /* this part only works if x is real and positive; only used with b t_INT */
     176          21 :     pari_sp av2 = avma;
     177          21 :     GEN q, v, c, s = real_1(l), t = real_0(l);
     178        2226 :     for (k = n-1; k >= 0; k--)
     179             :     { /* T = (a+k)*(a1+k) = a*a1 + k(a+a1) + k^2 = previous(T) - S - 2k + 1 */
     180        2226 :       GEN p1 = gdiv(T, mulss(-n, k+1));
     181        2226 :       s = gprec_wtrunc(gaddgs(gmul(p1,s), 1), l);
     182        2226 :       t = gprec_wtrunc(gadd(gmul(p1,t), gaddgs(a,k)), l);
     183        2226 :       if (!k) break;
     184        2205 :       T = gsubgs(gsub(T, S), 2*k-1);
     185        2205 :       if (gc_needed(av2,3)) gerepileall(av2, 3, &s,&t,&T);
     186             :     }
     187          21 :     q = utor(n, l);
     188          21 :     zf = gpow(utoi(n), gneg_i(a), l);
     189          21 :     u = gprec_wensure(gmul(zf, s), l);
     190          21 :     v = gprec_wensure(gmul(zf, gdivgs(t,-n)), l);
     191             :     for(;;)
     192         490 :     {
     193         511 :       GEN p1, e, f, d = real_1(l), qmb = gsub(q,b);
     194             :       pari_sp av3;
     195         511 :       c = divur(5,q); if (expo(c) >= -1) c = real2n(-1, l);
     196         511 :       p1 = subsr(1, divrr(x,q)); if (cmprr(c,p1) > 0) c = p1;
     197         511 :       togglesign(c); av3 = avma;
     198         511 :       e = u;
     199         511 :       f = v;
     200         511 :       for(k = 1;; k++)
     201       47126 :       {
     202       47637 :         GEN w = gadd(gmul(gaddgs(a,k-1),u), gmul(gaddgs(qmb,1-k),v));
     203       47637 :         u = gmul(divru(q,k),v);
     204       47637 :         v = gdivgu(w, k);
     205       47637 :         d = mulrr(d, c);
     206       47637 :         e = gadd(e, gmul(d,u));
     207       47637 :         f = gadd(f, p1 = gmul(d,v));
     208       47637 :         if (gequal0(p1) || gexpo(p1) - gexpo(f) <= 1-prec2nbits(precision(p1)))
     209             :           break;
     210       47126 :         if (gc_needed(av3,3)) gerepileall(av3,5,&u,&v,&d,&e,&f);
     211             :       }
     212         511 :       u = e;
     213         511 :       v = f;
     214         511 :       q = mulrr(q, addrs(c,1));
     215         511 :       if (expo(x) - expo(subrr(q,x)) >= bit) break;
     216         490 :       gerepileall(av2, 3, &u,&v,&q);
     217             :     }
     218             :   }
     219             :   else
     220             :   { /* this part works for large complex x */
     221          14 :     GEN zz = gneg_i(ginv(x)), s = gen_1;
     222          14 :     zf = gpow(x, gneg_i(a), l);
     223        1589 :     for (k = n-1; k >= 0; k--)
     224             :     {
     225        1589 :       s = gaddsg(1, gmul(gmul(T, gdivgu(zz,k+1)), s));
     226        1589 :       if (!k) break;
     227        1575 :       T = gsubgs(gsub(T, S), 2*k-1);
     228             :     }
     229          14 :     u = gmul(s, zf);
     230             :   }
     231          35 :   return gprec_wtrunc(u, prec);
     232             : }
     233             : GEN
     234          14 : hyperu(GEN a, GEN b, GEN x, long prec)
     235          14 : { pari_sp av = avma; return gerepilecopy(av, hyperu_i(a,b,x,prec)); }
     236             : 
     237             : static GEN
     238         364 : mkendpt(GEN z, GEN a)
     239             : {
     240         364 :   a = real_i(a);
     241         364 :   if (gcmpgs(a,-1) <= 0) pari_err_IMPL("hypergeom for these parameters");
     242         364 :   return (gcmpgs(a,1) >= 0 || gequal0(a))? z: mkvec2(z, a);
     243             : }
     244             : 
     245             : /*z != 0 */
     246             : static GEN
     247          56 : F20(GEN a, GEN b, GEN z, long prec)
     248             : {
     249             :   GEN U;
     250          56 :   z = gneg_i(z); U = hyperu_i(a, gadd(a, gsubsg(1, b)), ginv(z), prec);
     251          56 :   return gmul(U, gpow(z, gneg(a), prec));
     252             : }
     253             : 
     254             : static GEN F21(GEN a, GEN b, GEN c, GEN z, long prec);
     255             : 
     256             : static GEN
     257       17472 : fF31(void *E, GEN t)
     258             : {
     259       17472 :   pari_sp av = avma;
     260       17472 :   GEN a1 = gel(E,1), b = gel(E,2), c = gel(E,3), d = gel(E,4), z = gel(E,5);
     261       17472 :   long prec = precision(t);
     262       17472 :   return gerepileupto(av, gmul(gmul(gexp(gneg(t), prec), gpow(t, a1, prec)),
     263             :                                F21(b, c, d, gmul(t, z), prec)));
     264             : }
     265             : /* F31(a,b,c;d,z) = \int_0^oo\exp(-t)t^{a-1}F_{21}(b,c;d,tz) / gamma(a) */
     266             : static GEN
     267          56 : F31(GEN a, GEN b, GEN c, GEN d, GEN z, long prec)
     268             : {
     269             :   GEN p1, p2, a1;
     270          56 :   if (gcmp(real_i(a), real_i(b)) < 0) swap(a,b);
     271          56 :   if (gcmp(real_i(a), real_i(c)) < 0) swap(a,c);
     272          56 :   if (gsigne(real_i(a)) <= 0) pari_err_IMPL("F31 with a, b, and c <= 0");
     273          56 :   a1 = gsubgs(a,1);
     274          56 :   p1 = mkendpt(gen_0, a1);
     275          56 :   p2 = mkvec2(mkoo(), gen_1);
     276          56 :   return gdiv(intnum(mkvecn(5, a1, b, c, d, z), fF31, p1, p2, NULL, prec),
     277             :               ggamma(a, prec));
     278             : }
     279             : 
     280             : /* F32(a,b,c;d,e;z)=\int_0^1,x^{c-1}(1-x)^{e-c-1}F21(a,b,d,x*z)*gamma(e)/(gamma(c)*gamma(e-c)) */
     281             : static GEN
     282       11298 : fF32(void *E, GEN x)
     283             : {
     284       11298 :   pari_sp av = avma;
     285       11298 :   GEN c1 = gel(E,1), ec = gel(E,2), a = gel(E,3), b = gel(E,4);
     286       11298 :   GEN d = gel(E,5), z = gel(E,6), T;
     287       11298 :   long prec = precision(x);
     288       11298 :   T = F21(a, b, d, gmul(x, z), prec);
     289       11298 :   if (!gequal0(c1)) T = gmul(T, gpow(x,c1,prec));
     290       11298 :   if (!gequal0(ec)) T = gmul(T, gpow(gsubsg(1,x),ec,prec));
     291       11298 :   return gerepileupto(av,T);
     292             : }
     293             : 
     294             : static GEN
     295          42 : myint32(GEN a, GEN b, GEN c, GEN d, GEN e, GEN z, long prec)
     296             : {
     297          42 :   GEN c1 = gsubgs(c,1), c2 = gsub(e, gaddgs(c,1));
     298          42 :   GEN p0 = mkendpt(gen_0, c1);
     299          42 :   GEN p1 = mkendpt(gen_1, c2);
     300          42 :   return intnum(mkvecn(6, c1, c2, a,b,d, z), fF32, p0, p1, NULL, prec);
     301             : }
     302             : 
     303             : static void
     304          91 : check_hyp1(GEN x)
     305             : {
     306          91 :   if (gsigne(real_i(x)) <= 0)
     307           7 :     pari_err_DOMAIN("hypergeom","real(vecsum(D)-vecsum(N))", "<=", gen_0, x);
     308          84 : }
     309             : 
     310             : /* 0 < re(z) < e ? */
     311             : static int
     312          63 : ok_F32(GEN z, GEN e)
     313          63 : { GEN x = real_i(z); return gsigne(x) > 0 && gcmp(e, x) > 0; }
     314             : 
     315             : /* |z| very close to 1 but z != 1 */
     316             : static GEN
     317          49 : F32(GEN N, GEN D, GEN z, long prec)
     318             : {
     319             :   GEN tmp, a,b,c, d,e, re;
     320          49 :   a = gel(N,1); d = gel(D,1);
     321          49 :   b = gel(N,2); e = gel(D,2);
     322          49 :   c = gel(N,3);
     323          49 :   if (gcmp(real_i(e), real_i(d)) < 0) swap(e,d);
     324          49 :   re = real_i(e);
     325          49 :   if (!ok_F32(c, re))
     326             :   {
     327           7 :     if (ok_F32(b, re))  {swap(b,c);}
     328           7 :     else if (ok_F32(a, re)) {swap(a,c);}
     329           7 :     else pari_err_IMPL("3F2 for these arguments");
     330             :   }
     331          42 :   tmp = gdiv(ggamma(e, prec), mulgamma2(c, gsub(e,c), prec));
     332          42 :   return gmul(tmp, myint32(a, b, c, d, e, z, prec));
     333             : }
     334             : 
     335             : static GEN
     336      114023 : poch(GEN a, long n, long prec)
     337             : {
     338      114023 :   GEN S = real_1(prec);
     339             :   long j;
     340     1387288 :   for (j = 0; j < n; j++) S = gmul(S, gaddsg(j, a));
     341      114023 :   return S;
     342             : }
     343             : static GEN
     344         140 : vpoch(GEN a, long n)
     345             : {
     346         140 :   GEN v = cgetg(n+1, t_VEC);
     347             :   long j;
     348         140 :   gel(v,1) = a;
     349        1260 :   for (j = 1; j < n; j++) gel(v,j+1) = gmul(gel(v,j), gaddsg(j, a));
     350         140 :   return v;
     351             : }
     352             : static GEN
     353          98 : RgV_vpoch(GEN N, long n)
     354             : {
     355             :   long i, l;
     356          98 :   GEN v = cgetg_copy(N, &l);
     357         182 :   for (i = 1; i < l; i++) gel(v,i) = vpoch(gel(N,i), n);
     358          98 :   return v;
     359             : }
     360             : 
     361             : static GEN
     362      113771 : Npoch(GEN a, long n) { return gnorm(poch(a, n, LOWDEFAULTPREC)); }
     363             : static GEN
     364       38367 : Npochden(GEN a, long n)
     365             : {
     366       38367 :   GEN L = Npoch(a, n), r = ground(real_i(a));
     367       38367 :   if (signe(r) <= 0)
     368             :   {
     369        6237 :     GEN t = gnorm(gsub(a, r));
     370        6237 :     if (gcmpgs(t,1) > 0) L = gmul(L, t);
     371             :   }
     372       38367 :   return L;
     373             : }
     374             : 
     375             : /* |x + z|^2 */
     376             : static GEN
     377      111615 : normpol2(GEN z)
     378             : {
     379      111615 :   GEN a = real_i(z), b = imag_i(z);
     380      111615 :   return deg2pol_shallow(gen_1, gmul2n(a,1), gadd(gsqr(a),gsqr(b)), 0);
     381             : }
     382             : /* \prod |x + v[i]|^2 */
     383             : static GEN
     384       75166 : vnormpol2(GEN v)
     385             : {
     386       75166 :   long i, l = lg(v);
     387             :   GEN P;
     388       75166 :   if (l == 1) return pol_1(0);
     389       75166 :   P = normpol2(gel(v,1));
     390      111615 :   for (i = 2; i < l; i++) P = RgX_mul(P, normpol2(gel(v,i)));
     391       75166 :   return P;
     392             : }
     393             : 
     394             : static long
     395       37583 : precFtaylor(GEN N, GEN D, GEN z, long *pmi)
     396             : {
     397       37583 :   GEN v, ma, P = vnormpol2(D), Q = vnormpol2(N), Nz = gnorm(z);
     398       37583 :   double wma, logNz = (gexpo(Nz) < -27)? -27: dbllog2(Nz) / 2;
     399       37583 :   long pr = LOWDEFAULTPREC, prec = precision(z);
     400       37583 :   long lN = lg(N), lD = lg(D), mi, j, i, lv;
     401             : 
     402       37583 :   P = RgX_shift_shallow(P, 2);
     403             :   /* avoid almost cancellation of leading coeff if |z| ~ 1 */
     404       37583 :   if (!prec || fabs(logNz) > 1e-38) Q = RgX_Rg_mul(Q,Nz);
     405      111587 :   for (j = 1, ma = NULL; j < lN; ++j)
     406             :   {
     407       74004 :     GEN Nj = gel(N,j);
     408       74004 :     if (isint(Nj,&Nj) && signe(Nj) <= 0
     409        5215 :         && (!ma || abscmpii(ma, Nj) < 0)) ma = Nj;
     410             :   }
     411             :   /* use less sharp fujiwara_bound_real(,1) ? */
     412       37583 :   v = ground(realroots(gsub(P,Q), mkvec2(gen_0,mkoo()), pr));
     413       37583 :   v = ZV_to_zv(v); lv = lg(v);
     414       37583 :   if (ma)
     415             :   {
     416        5047 :     long sma = is_bigint(ma)? LONG_MAX: maxss(labs(itos(ma))-1, 1);
     417       10738 :     for (i = 1; i < lv; ++i) v[i] = maxss(minss(sma, v[i]), 1);
     418             :   }
     419       75922 :   for (i = 1, wma = 0., mi = 0; i < lv; ++i)
     420             :   {
     421       38339 :     GEN t1 = gen_1, t2 = gen_1;
     422       38339 :     long u = v[i];
     423             :     double t;
     424       38339 :     mi = maxss(mi, u);
     425      113743 :     for (j = 1; j < lN; j++) t1 = gmul(t1, Npoch(gel(N,j), u));
     426       76706 :     for (j = 1; j < lD; j++) t2 = gmul(t2, Npochden(gel(D,j), u));
     427       38339 :     t = dbllog2(gdiv(t1,t2)) / 2 + u * logNz - dbllog2(mpfactr(u,pr));
     428       38339 :     wma = maxdd(wma, t); /* t ~ log2 | N(u)/D(u) z^u/u! | */
     429             :   }
     430             :   /* make up for exponential decrease in exp() */
     431       37583 :   if (gsigne(real_i(z)) < 0) wma -= gtodouble(real_i(z)) / M_LN2;
     432       37583 :   *pmi = mi; return ceil(wma/BITS_IN_LONG) + 1;
     433             : }
     434             : 
     435             : static GEN
     436       24374 : Ftaylor(GEN N, GEN D, GEN z, long prec)
     437             : {
     438             :   pari_sp av;
     439             :   GEN C, S;
     440       24374 :   long i, j, ct, lN = lg(N), lD = lg(D), pradd, mi, bitmin, tol;
     441       24374 :   pradd = precFtaylor(N, D, z, &mi);
     442       24374 :   if (pradd > 0)
     443             :   {
     444       24374 :     prec += pradd;
     445       24374 :     N = gprec_wensure(N, prec);
     446       24374 :     D = gprec_wensure(D, prec);
     447       24374 :     z = gprec_wensure(z, prec);
     448             :   }
     449       24374 :   bitmin = -(prec2nbits(prec) + 10);
     450       24374 :   S = C = real_1(prec); ct = 0; j = 0; tol = 0;
     451       24374 :   av = avma;
     452             :   for(;;)
     453     7238797 :   {
     454     7263171 :     GEN a = gen_1, b = gen_1;
     455    21117545 :     for (i = 1; i < lN; i++) a = gmul(a, gaddsg(j, gel(N,i)));
     456    14527325 :     for (i = 1; i < lD; i++) b = gmul(b, gaddsg(j, gel(D,i)));
     457     7263171 :     C = gmul(C, gmul(gdiv(a, gmulsg(j+1, b)), z));
     458     7263171 :     if (gequal0(C)) break;
     459     7262891 :     if (j > mi) tol = gequal0(S)? 0: gexpo(C) - gexpo(S);
     460     7262891 :     S = gadd(S, C); ++j;
     461     7262891 :     if (j > mi)
     462     7052114 :     { if (tol > bitmin) ct = 0; else if (++ct >= lN+lD-2) break; }
     463     7238797 :     if (gc_needed(av, 1)) gerepileall(av, 2, &S, &C);
     464             :   }
     465       24374 :   return S;
     466             : }
     467             : 
     468             : static GEN
     469        4704 : bind(GEN a, GEN b, GEN c, long ind)
     470             : {
     471        4704 :   switch(ind)
     472             :   {
     473         112 :     case 1: case 2: return gsub(c, b);
     474          56 :     case 5: case 6: return gsub(gaddsg(1, a), c);
     475        4536 :     default: return b;
     476             :   }
     477             : }
     478             : static GEN
     479      122024 : cind(GEN a, GEN b, GEN c, long ind)
     480             : {
     481      122024 :   switch(ind)
     482             :   {
     483       58716 :     case 1: case 6: return gaddsg(1, gsub(a,b));
     484       58716 :     case 4: case 5: return gsub(gaddsg(1, gadd(a,b)), c);
     485        4592 :     default: return c;
     486             :   }
     487             : }
     488             : static GEN
     489      115640 : zind(GEN z, long ind)
     490             : {
     491      115640 :   switch(ind)
     492             :   {
     493       24738 :     case 1: return ginv(gsubsg(1, z));
     494       29414 :     case 2: return gdiv(z, gsubgs(z, 1));
     495       16107 :     case 4: return gsubsg(1, z);
     496       16107 :     case 5: return gsubsg(1, ginv(z));
     497       24738 :     case 6: return ginv(z);
     498             :   }
     499        4536 :   return z;
     500             : }
     501             : 
     502             : /* z not 0 or 1, c not a nonpositive integer */
     503             : static long
     504       29330 : F21ind(GEN a, GEN b, GEN c, GEN z, long bit)
     505             : {
     506       29330 :   GEN v = const_vec(6, mkoo());
     507       29330 :   long ind = 0;
     508       29330 :   const long LD = LOWDEFAULTPREC;
     509       29330 :   if (!isnegint(cind(a,b,c, 1))) gel(v,1) = gabs(zind(z,1), LD);
     510       29330 :   gel(v,2) = gabs(zind(z,2), LD);
     511       29330 :   gel(v,3) = gabs(z, LD);
     512       29330 :   if (!isnegint(cind(a,b,c, 4))) gel(v,4) = gabs(zind(z,4), LD);
     513       29330 :   if (!isnegint(cind(a,b,c, 5))) gel(v,5) = gabs(zind(z,5), LD);
     514       29330 :   if (!isnegint(cind(a,b,c, 6))) gel(v,6) = gabs(zind(z,6), LD);
     515       29330 :   ind = vecindexmin(v); /* |znew| <= 1; close to 1 ? */
     516       29330 :   return (gexpo(gsubgs(gel(v,ind),1)) > -maxss(bit / 4, 32))? -ind: ind;
     517             : }
     518             : static GEN
     519       13258 : mul4(GEN a, GEN b, GEN c, GEN d) { return gmul(a,gmul(b, gmul(c, d))); }
     520             : static GEN
     521       56154 : mul3(GEN a, GEN b, GEN c) { return gmul(a,gmul(b, c)); }
     522             : 
     523             : /* (1 - zt)^a t^b (1-t)^c */
     524             : static GEN
     525       56042 : fF212(void *E, GEN t)
     526             : {
     527       56042 :   GEN z = gel(E,1), a = gel(E,2), b = gel(E,3), c = gel(E,4);
     528       56042 :   GEN u = gsubsg(1, gmul(z, t));
     529       56042 :   long prec = precision(t);
     530       56042 :   return mul3(gpow(u, a, prec), gpow(t, b, prec), gpow(gsubsg(1,t), c, prec));
     531             : }
     532             : 
     533             : /* (1 - zt)^a T(1-zt) t^b (1-t)^c */
     534             : static GEN
     535       13258 : fF21neg2(void *E, GEN t)
     536             : {
     537       13258 :   GEN z = gel(E,1), a = gel(E,2), b = gel(E,3), c = gel(E,4), T = gel(E,5);
     538       13258 :   GEN u = gsubsg(1, gmul(z, t));
     539       13258 :   long prec = precision(t);
     540       13258 :   return mul4(gsubst(T, 0, u), gpow(u, a, prec), gpow(t, b, prec),
     541             :               gpow(gsubsg(1,t), c, prec));
     542             : }
     543             : 
     544             : /* N >= 1 */
     545             : static GEN
     546          28 : F21lam(long N, GEN a, GEN c)
     547             : {
     548             :   long i;
     549          28 :   GEN C = vecbinomial(N), S = cgetg(N+2, t_VEC);
     550          28 :   GEN vb = vpoch(gsub(c,a), N), va = vpoch(a, N);
     551          28 :   gel(S,1) = gel(va,N);
     552          28 :   for (i = 1; i < N; i++) gel(S,i+1) = mul3(gel(C,i+1), gel(vb,i), gel(va,N-i));
     553          28 :   gel(S,i+1) = gel(vb,N); return RgV_to_RgX(S,0);
     554             : }
     555             : 
     556             : /* F(-m,b; c; z), m >= 0 */
     557             : static GEN
     558        4732 : F21finitetaylor(long m, GEN b, GEN c, GEN z, long prec)
     559             : {
     560             :   pari_sp av;
     561             :   GEN C, S;
     562             :   long j, ct, pradd, mi, bitmin, mb;
     563        4732 :   if (isnegint2(b, &mb) && mb < m) { b = utoineg(m); m = mb; }
     564        4732 :   pradd = precFtaylor(mkvec2(utoineg(m), b), mkvec(c), z, &mi);
     565        4732 :   if (pradd > 0)
     566             :   {
     567        4732 :     prec += pradd;
     568        4732 :     b = gprec_wensure(b, prec);
     569        4732 :     c = gprec_wensure(c, prec);
     570        4732 :     z = gprec_wensure(z, prec);
     571             :   }
     572        4732 :   bitmin = -(prec2nbits(prec) + 10);
     573        4732 :   C = real_1(prec); S = C; ct = 0;
     574        4732 :   av = avma;
     575      133168 :   for(j = 0; j < m; ++j)
     576             :   {
     577      128436 :     C = gmul(C, gdiv(gmulsg(j-m, gaddsg(j, b)), gmulsg(j+1, gaddsg(j, c))));
     578      128436 :     C = gmul(C, z);
     579      128436 :     if (j > mi && !gequal0(S))
     580        7084 :     { if (gexpo(C) - gexpo(S) > bitmin) ct = 0; else if (++ct == 3) break; }
     581      128436 :     S = gadd(S, C);
     582      128436 :     if (gc_needed(av, 1)) gerepileall(av, 2, &S, &C);
     583             :   }
     584        4732 :   return S;
     585             : }
     586             : 
     587             : /* c not a nonpositive integer */
     588             : static GEN
     589         112 : F21finite_i(long m, GEN b, GEN c, GEN z, GEN B, GEN C, GEN coe, long prec)
     590             : {
     591         112 :   return mul3(poch(B, m, prec), gdiv(gpowgs(coe, m), poch(C, m, prec)),
     592             :               F21finitetaylor(m, b, c, z, prec));
     593             : }
     594             : 
     595             : /* F(-m,b; c; z), m >= 0; c not a nonpositive integer */
     596             : static GEN
     597        4732 : F21finite(long m, GEN b, GEN c, GEN z, long prec)
     598             : {
     599        4732 :   GEN a = stoi(-m), b1 = b, c1 = c, z1;
     600        4732 :   long ind = F21ind(a, b, c, z, prec2nbits(prec)), inda = labs(ind);
     601        4732 :   z1 = zind(z, inda);
     602        4732 :   if (ind < 0)
     603             :   {
     604        4704 :     b1 = bind(a, b, c, inda);
     605        4704 :     c1 = cind(a, b, c, inda); /* not a nonpositive integer */
     606             :   }
     607        4732 :   switch (inda)
     608             :   {
     609          28 :     case 1: return F21finite_i(m, b1, c1, z1, b, c, gsubsg(1,z), prec);
     610          84 :     case 2: return gmul(gpowgs(gsubsg(1,z), m),
     611             :                         F21finitetaylor(m, b1, c, z1, prec));
     612          28 :     case 4: return F21finite_i(m, b1, c1, z1, gsub(c,b), c, gen_1, prec);
     613          28 :     case 5: return F21finite_i(m, b1, c1, z1, gsub(c,b), c, z, prec);
     614          28 :     case 6: return F21finite_i(m, b1, c1, z1, b, c, gneg(z), prec);
     615        4536 :     default:return F21finitetaylor(m, b1, c1, z, prec);
     616             :   }
     617             : }
     618             : 
     619             : /**********************************************************/
     620             : 
     621             : static GEN
     622         140 : multgam(GEN a, GEN b, GEN c, GEN d, long prec)
     623             : {
     624         140 :   if (isnegint(c) || isnegint(d)) return gen_0;
     625         140 :   return gdiv(mulgamma2(a, b, prec), mulgamma2(c, d, prec));
     626             : }
     627             : 
     628             : static GEN
     629         112 : intnumsplit(void *E, GEN (*f)(void*, GEN), GEN a, GEN b, GEN z, long prec)
     630             : {
     631         112 :   if (!z) return intnum(E, f, a, b, NULL, prec);
     632           0 :   return gadd(intnum(E, f, a, z, NULL, prec),
     633             :               intnum(E, f, z, b, NULL, prec));
     634             : }
     635             : /* z != 1 */
     636             : static GEN
     637         112 : myint21(void *E, GEN (*f)(void*, GEN), long prec)
     638             : {
     639         112 :   GEN z = gel(E,1), a = real_i(gel(E,2)), b = gel(E,3), c = gel(E,4);
     640         112 :   GEN pz = NULL, p0 = mkendpt(gen_0, b), p1 = mkendpt(gen_1, c);
     641         112 :   if (gcmpgs(a, 1) <= 0 && is0(imag_i(z), 10))
     642             :   {
     643             :     GEN r;
     644           0 :     pz = ginv(z); r = real_i(pz);
     645           0 :     if (gsigne(r) <= 0 || gcmp(r, gen_1) >= 0) pz = NULL;
     646             :   }
     647         112 :   if (pz) pz = mkendpt(pz,a);
     648         112 :   else if (gcmpgs(a,-1) <= 0) prec += ((gexpo(a)+1)>>1) * EXTRAPREC64;
     649         112 :   return intnumsplit(E, f, p0, p1, pz, prec);
     650             : }
     651             : 
     652             : /* Algorithm used for F21(a,b;c;z)
     653             : Basic transforms:
     654             :   1: (c-b,1+a-b,1/(1-z))
     655             :   2: (c-b,c,z/(z-1))
     656             :   3: (b,c,z)
     657             :   4: (b,b-c+a+1,1-z)
     658             :   5: (1+a-c,b-c+a+1,1-1/z)
     659             :   6: (1+a-c,1+a-b,1/z)
     660             : 
     661             : F21: calls F21_i and increase prec if too much cancellation
     662             : F21_i: c is not a non-positive integer
     663             : - z ~ 0 or 1: return special value
     664             : - if a, b, c-b or c-a a non-positive integer: use F21finite
     665             : - compute index, value of z
     666             :    if |z| < 1-epsilon return F21taylorind
     667             :    if Re(b)<=0, swap and/or recurse
     668             :    so may assume Re(b)>0 and Re(a)>=Re(b) and integrate.
     669             : 
     670             : F21finite:
     671             : - compute index, value of z
     672             : - call F21finitetaylor
     673             : 
     674             : F21ind: find best index (1 to 6, -1 to -6 if |z| < 1-epsilon)
     675             : F21finitetaylor: a or b in Z_{<=0}; calls precFtaylor
     676             : 
     677             : F21taylorind: in case 2, may lose accuracy, possible bug.
     678             : - calls F21taylor[1456] or F21taylor
     679             : 
     680             : F21taylor: calls Ftaylor / gamma: may lose accuracy
     681             : 
     682             : FBaux1: F21taylor twice
     683             : FBaux2: F21finitelim + F21taylorlim
     684             : F21taylor[45]: if c-(a+b) integer, calls FBaux1, else calls FBaux2
     685             : F21taylor[16]: if b-a integer, calls FBaux1, else calls FBaux2
     686             : 
     687             : F21taylorlim: calls precFtaylor then compute
     688             : F21finitelim: direct */
     689             : static GEN F21taylorind(GEN a, GEN b, GEN c, GEN z, long ind, long prec);
     690             : /* c not a nonpositive integer */
     691             : static GEN
     692       30219 : F21_i(GEN a, GEN b, GEN c, GEN z, long prec)
     693             : {
     694             :   GEN res;
     695       30219 :   long m, ind, prec2, bitprec = prec2nbits(prec);
     696       30219 :   if (is0(imag_i(z), bitprec)) z = real_i(z);
     697       30219 :   if (is0(z, bitprec)) return real_1(prec);
     698       29365 :   if (gequal1(z))
     699             :   {
     700          35 :     GEN x = gsub(c, gadd(a, b)); check_hyp1(x);
     701          28 :     return multgam(c, x, gsub(c,a), gsub(c,b), prec);
     702             :   }
     703       29330 :   if (isnegint2(b, &m)) return F21finite(m, a, c, z, prec);
     704       29134 :   if (isnegint2(a, &m)) return F21finite(m, b, c, z, prec);
     705       24682 :   if (isnegint(gsub(c, b))) swap(a, b);
     706       24682 :   if (isnegint2(gsub(c, a), &m))
     707             :   {
     708          84 :     GEN x = gpow(gsubsg(1, z), gneg(gaddsg(m, b)), prec);
     709          84 :     return gmul(x, F21finite(m, gsub(c, b), c, z, prec));
     710             :   }
     711             :   /* Here a, b, c, c-a, c-b are not nonpositive integers */
     712       24598 :   ind = F21ind(a, b, c, z, bitprec);
     713       24598 :   prec2 = prec + EXTRAPREC64;
     714       24598 :   a = gprec_wensure(a,prec2);
     715       24598 :   b = gprec_wensure(b,prec2);
     716       24598 :   c = gprec_wensure(c,prec2);
     717       24598 :   z = gprec_wensure(z,prec2);
     718       24598 :   if (ind < 0) return gmul(ggamma(c, prec), F21taylorind(a,b,c, z, ind, prec));
     719         112 :   if (gsigne(real_i(b)) <= 0)
     720             :   {
     721          28 :     if (gsigne(real_i(a)) <= 0)
     722             :     {
     723             :       GEN p1,p2;
     724           0 :       if (gcmp(real_i(b), real_i(a)) < 0) swap(a,b);
     725             :       /* FIXME: solve recursion as below with F21auxpol */
     726           0 :       p1 = gmul(gsubsg(1, z), F21_i(a, gaddsg(1,b), c, z, prec));
     727           0 :       p2 = gmul(gmul(gsubsg(1, gdiv(a,c)), z),
     728             :                 F21_i(a, gaddsg(1,b), gaddsg(1,c), z, prec));
     729           0 :       return gadd(p1, p2);
     730             :     }
     731          28 :     swap(a,b);
     732             :   }
     733         112 :   if (gcmp(real_i(a), real_i(b)) < 0 && gsigne(real_i(a)) > 0) swap(a,b);
     734             :   /* Here real(b) > 0 and either real(a) <= 0 or real(a) > real(b) */
     735         112 :   if (gcmp(real_i(c), real_i(b)) <= 0)
     736             :   {
     737          28 :     long N = 1 + itos(gfloor(gsub(real_i(b),real_i(c)))); /* >= 1 */
     738          28 :     GEN T = F21lam(N, a, c), c0 = c;
     739             :     void *E;
     740          28 :     c = gaddsg(N,c);
     741          28 :     E = (void*)mkvec5(z, gsubsg(-N,a), gsubgs(b,1), gsubgs(gsub(c,b),1), T);
     742          28 :     res = gdiv(myint21(E, fF21neg2, prec2), poch(c0, N, prec));
     743             :   }
     744             :   else
     745             :   {
     746          84 :     void *E = (void*)mkvec4(z, gneg(a), gsubgs(b,1), gsubgs(gsub(c,b),1));
     747          84 :     res = myint21(E, fF212, prec2);
     748             :   }
     749         112 :   return gmul(multgam(gen_1, c, b, gsub(c,b), prec), res);
     750             : }
     751             : 
     752             : /* c not a non-positive integer */
     753             : static GEN
     754       30163 : F21(GEN a, GEN b, GEN c, GEN z, long prec)
     755             : {
     756       30163 :   GEN res = F21_i(a, b, c, z, prec);
     757       30156 :   long ex = labs(gexpo(res)), bitprec = prec2nbits(prec);
     758       30156 :   if (ex > bitprec)
     759             :   {
     760          56 :     prec = nbits2prec(ex + bitprec);
     761          56 :     res = F21_i(gprec_wensure(a,prec), gprec_wensure(b,prec),
     762             :                 gprec_wensure(c,prec), gprec_wensure(z,prec), prec);
     763             :   }
     764       30156 :   return res;
     765             : }
     766             : 
     767             : static GEN
     768       23170 : F21taylor(GEN a, GEN b, GEN c, GEN z, long prec)
     769             : {
     770       23170 :   pari_sp av = avma;
     771       23170 :   GEN r = gdiv(Ftaylor(mkvec2(a,b), mkvec(c), z, prec), ggamma(c, prec));
     772       23170 :   return gerepileupto(av, r);
     773             : }
     774             : 
     775             : static GEN
     776        8477 : F21taylorlim(GEN N, long m, GEN z, GEN Z, long ind, long prec)
     777             : {
     778             :   pari_sp av;
     779             :   GEN C, P, S, a, b;
     780        8477 :   long j, ct, pradd, mi, fl, bitmin, tol, si = (ind == 5 || ind == 6)? -1: 1;
     781        8477 :   pradd = precFtaylor(N, mkvec(stoi(m + 1)), z, &mi);
     782        8477 :   if (pradd)
     783             :   {
     784        8477 :     prec += pradd;
     785        8477 :     N = gprec_wensure(N, prec);
     786        8477 :     z = gprec_wensure(z, prec);
     787        8477 :     Z = gprec_wensure(Z, prec);
     788             :   }
     789        8477 :   av = avma; a = gel(N,1); b = gel(N,2);
     790        8477 :   bitmin = -(prec2nbits(prec) + 10);
     791        8477 :   P = glog(Z, prec); if (ind == 4 || ind == 5) P = gneg(P);
     792        8477 :   P = gadd(P, gsub(gpsi(stoi(m+1), prec), mpeuler(prec)));
     793        8477 :   P = gsub(P, gadd(gpsi(a, prec), gpsi(si == -1? gsubsg(1, b): b, prec)));
     794        8477 :   C = real_1(prec); ct = 0; tol = 0; S = P;
     795        8477 :   for(j = 0, fl = 1;;)
     796      993465 :   {
     797     1001942 :     GEN v1 = gaddgs(a, j), v2 = gaddgs(b, j);
     798     1001942 :     long J = (j+1) * (j+1+m);
     799     1001942 :     C = gmul(C, gdivgs(gmul(z, v1), J));
     800     1001942 :     if (gequal0(v2)) fl = 0; else C = gmul(C, v2);
     801     1001942 :     if (j > mi) tol = gequal0(S) ? 0 : gexpo(C) - gexpo(S);
     802     1001942 :     if (fl)
     803             :     {
     804     1000249 :       P = gadd(P, gsub(uutoQ(2*j+2+m, J), gadd(ginv(v1), ginv(v2))));
     805     1000249 :       S = gadd(S, gmul(C, P));
     806             :     }
     807             :     else
     808        1693 :       S = (si == 1)? gadd(S, C): gsub(S, C);
     809     1001942 :     if (++j > mi) { if (tol > bitmin) ct = 0; else if (++ct == 3) break; }
     810      993465 :     if (gc_needed(av, 1)) gerepileall(av, 3, &S, &C, &P);
     811             :   }
     812        8477 :   return gdiv(S, mpfact(m));
     813             : }
     814             : 
     815             : /* N = [a,b]; (m-1)! sum_{1 <= k < m} (-z)^k (a)_k (b)_k / k! (m-k)!*/
     816             : static GEN
     817        8477 : F21finitelim(GEN N, long m, GEN z, long prec)
     818             : {
     819             :   GEN C, S, a, b;
     820             :   long j;
     821        8477 :   if (!m) return gen_0;
     822          91 :   a = gel(N,1);
     823          91 :   b = gel(N,2);
     824          91 :   S = C = real_1(prec + EXTRAPRECWORD);
     825         147 :   for (j = 1; j < m; j++)
     826             :   {
     827          56 :     GEN v1 = gaddsg(j-1, a), v2 = gaddsg(j-1, b);
     828          56 :     C = gdivgs(gmul(C, gmul(gmul(v1, v2), z)), j*(j-m));
     829          56 :     S = gadd(S, C);
     830             :   }
     831          91 :   return gmul(S, mpfact(m-1));
     832             : }
     833             : 
     834             : static GEN
     835        5796 : OK_gadd(GEN x, GEN y, long prec0, long *pprec,
     836             :         GEN *d1,GEN *d2,GEN *d3,GEN *d4,GEN *d5,GEN *d6,GEN *d7,GEN *d8)
     837             : {
     838        5796 :   long prec = *pprec;
     839        5796 :   GEN z = gadd(x,y);
     840        5796 :   if (!gequal0(z) && gexpo(z)-gexpo(x) >= prec2nbits(prec0) - prec2nbits(prec))
     841        4431 :     return z;
     842        1365 :   *pprec = prec = precdbl(prec);
     843        1365 :   *d1 = gprec_wensure(*d1, prec); *d2 = gprec_wensure(*d2, prec);
     844        1365 :   *d3 = gprec_wensure(*d3, prec); *d4 = gprec_wensure(*d4, prec);
     845        1365 :   *d5 = gprec_wensure(*d5, prec); *d6 = gprec_wensure(*d6, prec);
     846        1365 :   *d7 = gprec_wensure(*d7, prec); *d8 = gprec_wensure(*d8, prec);
     847        1365 :   return NULL;
     848             : }
     849             : 
     850             : static GEN
     851        4431 : FBaux1(GEN v1, GEN g1, GEN c1, GEN v2, GEN g2, GEN c2, GEN z, GEN bma,
     852             :        long prec0, long prec)
     853             : {
     854        4431 :   GEN pi = mppi(prec);
     855             :   for (;;)
     856        1365 :   {
     857        5796 :     GEN t1 = gdiv(c1, mulgammav2(g1, prec)), r1;
     858        5796 :     GEN t2 = gdiv(c2, mulgammav2(g2, prec)), r2, F;
     859        5796 :     r1 = gmul(t1, F21taylor(gel(v1,1), gel(v1,2), gel(v1,3), z, prec));
     860        5796 :     r2 = gmul(t2, F21taylor(gel(v2,1), gel(v2,2), gel(v2,3), z, prec));
     861        5796 :     F = OK_gadd(r1, r2, prec0, &prec, &c1,&c2, &g1,&g2, &v1,&v2, &z,&bma);
     862        5796 :     if (F) return gmul(F, gdiv(pi, gsin(gmul(pi, bma), prec)));
     863             :   }
     864             : }
     865             : 
     866             : static GEN
     867        8477 : FBaux2(GEN v1, GEN g1, GEN c1, long m, GEN z1, GEN c2, GEN g2, GEN v2, GEN z2,
     868             :        GEN Z2, long ind, long prec)
     869             : {
     870        8477 :   GEN t1 = gdiv(c1, mulgammav2(g1, prec)), r1;
     871        8477 :   GEN t2 = gdiv(c2, mulgammav2(g2, prec)), r2;
     872        8477 :   r1 = gmul(t1, F21finitelim(v1, m, z1, prec));
     873        8477 :   r2 = gmul(t2, F21taylorlim(v2, m, z2, Z2, ind, prec));
     874        8477 :   return gadd(r1, r2);
     875             : }
     876             : 
     877             : /* 1 / (1-z) */
     878             : static GEN
     879       12572 : F21taylor1(GEN a, GEN b, GEN c, GEN z, long prec)
     880             : {
     881       12572 :   GEN bma = gsub(b, a), coe1, coe2, z1, g1, g2, v1, v2, Z;
     882             :   long m;
     883       12572 :   if (!islong(bma,&m, prec))
     884             :   {
     885             :     GEN b1, c1, e1, c2;
     886        4207 :     b1 = gsub(c, b);
     887        4207 :     c1 = gsubsg(1, bma);
     888        4207 :     e1 = gsub(c, a);
     889        4207 :     coe1 = gpow(gsubsg(1, z), gneg(a), prec);
     890        4207 :     c2 = gaddgs(bma, 1);
     891        4207 :     coe2 = gneg(gpow(gsubsg(1, z), gneg(b), prec));
     892        4207 :     z1 = ginv(gsubsg(1, z));
     893        4207 :     return FBaux1(mkvec3(a,b1,c1), mkvec2(b,e1), coe1, mkvec3(b,e1,c2),
     894             :                   mkvec2(a,b1), coe2, z1, bma, prec, prec);
     895             :   }
     896        8365 :   if (m < 0) { swap(a,b); m = -m; }
     897        8365 :   Z = gsubsg(1, z);
     898        8365 :   coe1 = gpow(Z, gneg(a), prec);
     899        8365 :   v2 = g1 = mkvec2(gaddgs(a, m), gsub(c, a));
     900        8365 :   v1 = mkvec2(a, gsub(c, gaddgs(a, m)));
     901        8365 :   z1 = ginv(Z);
     902        8365 :   coe2 = gmul(coe1, gpowgs(z1, m)); if (m & 1L) coe2 = gneg(coe2);
     903        8365 :   g2 = mkvec2(a, gsub(c, gaddgs(a, m))); /* 15.8.9 */
     904        8365 :   return FBaux2(v1, g1, coe1, m, z1, coe2, g2, v2, z1, Z, 1, prec);
     905             : }
     906             : 
     907             : /* 1 - z */
     908             : static GEN
     909         210 : F21taylor4(GEN a, GEN b, GEN c, GEN z, long prec)
     910             : {
     911         210 :   GEN bma = gsub(c, gadd(a, b)), coe2, z1, z2, g1, g2, v1, v2;
     912             :   long m;
     913         210 :   if (!islong(bma,&m, prec))
     914             :   {
     915             :     GEN c1, a2, b2, c2;
     916         119 :     c1 = gsubsg(1, bma);
     917         119 :     a2 = gsub(c, a);
     918         119 :     b2 = gsub(c, b);
     919         119 :     c2 = gaddsg(1, bma);
     920         119 :     z1 = gsubsg(1, z); coe2 = gneg(gpow(z1, bma, prec));
     921         119 :     return FBaux1(mkvec3(a,b,c1), mkvec2(a2,b2), gen_1, mkvec3(a2,b2,c2),
     922             :                   mkvec2(a,b), coe2, z1, bma, prec, prec);
     923             :   }
     924          91 :   if (m < 0)
     925             :   {
     926          42 :     GEN F = F21taylor4(gaddgs(a,m), gaddgs(b,m), gaddgs(gadd(a,b), m), z, prec);
     927          42 :     return gmul(gpowgs(gsubsg(1,z), m), F);
     928             :   }
     929          49 :   v2 = g1 = mkvec2(gaddgs(a,m), gaddgs(b,m));
     930          49 :   v1 = g2 = mkvec2(a, b);
     931          49 :   z1 = gsubgs(z, 1);
     932          49 :   z2 = gneg(z1); coe2 = gpowgs(z1, m); /* 15.8.10 */
     933          49 :   return FBaux2(v1, g1, gen_1, m, z1, coe2, g2, v2, z2, z2, 4, prec);
     934             : }
     935             : 
     936             : /* 1 - 1/z */
     937             : static GEN
     938         112 : F21taylor5(GEN a, GEN b, GEN c, GEN z, long prec)
     939             : {
     940         112 :   GEN bma = gsub(c, gadd(a,b)), tmp, coe1, coe2, z1, g1, g2, v1, v2, z2;
     941             :   long m;
     942         112 :   if (!islong(bma,&m, prec))
     943             :   {
     944             :     GEN b1, c1, d1, e1, b2, c2;
     945          49 :     d1 = gsub(c, a);
     946          49 :     b1 = gsubsg(1, d1);
     947          49 :     c1 = gsubsg(1, bma);
     948          49 :     e1 = gsub(c, b);
     949          49 :     b2 = gsubsg(1, a);
     950          49 :     c2 = gaddsg(1, bma);
     951          49 :     coe1 = gpow(z, gneg(a), prec);
     952          49 :     coe2 = gneg(gmul(gpow(gsubsg(1, z), bma, prec), gpow(z, gneg(d1), prec)));
     953          49 :     z1 = gsubsg(1, ginv(z));
     954          49 :     return FBaux1(mkvec3(a,b1,c1), mkvec2(d1,e1), coe1,
     955             :                   mkvec3(d1,b2,c2), mkvec2(a, b), coe2, z1, bma, prec, prec);
     956             :   }
     957             :   /* c - (a + b) ~ m */
     958          63 :   if (m < 0)
     959             :   {
     960          28 :     tmp = F21taylor5(gaddgs(a,m), gaddgs(b,m), c, z, prec);
     961          28 :     return gmul(gpowgs(gsubsg(1, z), m), tmp);
     962             :   }
     963          35 :   g1 = mkvec2(gaddgs(a,m), gaddgs(b,m));
     964          35 :   v1 = mkvec2(a, gsubsg(1-m, b));
     965          35 :   v2 = mkvec2(gaddgs(a,m), gsubsg(1,b));
     966          35 :   z1 = gsubgs(ginv(z), 1);
     967          35 :   z2 = gneg(z1);
     968          35 :   g2 = mkvec2(a, b);
     969          35 :   coe1 = gpow(z, gneg(a), prec);
     970          35 :   coe2 = gmul(coe1, gpowgs(z2, m)); /* 15.8.11 */
     971          35 :   return FBaux2(v1, g1, coe1, m, z1, coe2, g2, v2, z2, z1, 5, prec);
     972             : }
     973             : 
     974             : /* 1 / z */
     975             : static GEN
     976          84 : F21taylor6(GEN a, GEN b, GEN c, GEN z, long prec)
     977             : {
     978          84 :   GEN bma = gsub(b, a), cma, am, coe1, coe2, z1, g1, g2, v1, v2, z2, Z;
     979             :   long m;
     980          84 :   if (!islong(bma,&m, prec))
     981             :   {
     982             :     GEN e1, e2, b1, b2, c1, c2;
     983          56 :     b1 = gaddgs(gsub(a,c), 1);
     984          56 :     c1 = gsubsg(1, bma);
     985          56 :     e1 = gsub(c,a);
     986          56 :     b2 = gaddgs(gsub(b,c), 1);
     987          56 :     c2 = gaddgs(bma, 1);
     988          56 :     e2 = gsub(c, b);
     989          56 :     coe1 = gpow(gneg(z), gneg(a), prec);
     990          56 :     coe2 = gneg(gpow(gneg(z), gneg(b), prec));
     991          56 :     z1 = ginv(z);
     992          56 :     return FBaux1(mkvec3(a,b1,c1), mkvec2(b,e1), coe1,
     993             :                   mkvec3(b,b2,c2), mkvec2(a,e2), coe2, z1, bma, prec, prec);
     994             :   }
     995             :   /* b - a ~ m */
     996          28 :   if (m < 0) { swap(a,b); m = -m; }
     997          28 :   cma = gsub(c, a); am = gaddgs(a,m);
     998          28 :   Z = gneg(z);
     999          28 :   coe1 = gpow(Z, gneg(a), prec);
    1000          28 :   coe2 = gdiv(coe1, gpowgs(z, m));
    1001          28 :   g1 = mkvec2(am, cma);
    1002          28 :   v1 = mkvec2(a, gsubsg(1, cma));
    1003          28 :   g2 = mkvec2(a, gsubgs(cma, m));
    1004          28 :   v2 = mkvec2(am, gsubsg(m+1, cma));
    1005          28 :   z2 = ginv(z);
    1006          28 :   z1 = gneg(z2); /* 15.8.8 */
    1007          28 :   return FBaux2(v1, g1, coe1, m, z1, coe2, g2, v2, z2, Z, 6, prec);
    1008             : }
    1009             : 
    1010             : /* (new b, new c, new z): given by bind, cind, zind
    1011             :  * case 1: (c-b,1+a-b,1/(1-z))
    1012             :  * case 2: (c-b,c,z/(z-1))
    1013             :  * case 3: (b,c,z)
    1014             :  * case 4: (b,b-c+a+1,1-z)
    1015             :  * case 5: (1+a-c,b-c+a+1,1-1/z)
    1016             :  * case 6: (1+a-c,1+a-b,1/z) */
    1017             : static GEN
    1018       24486 : F21taylorind(GEN a, GEN b, GEN c, GEN z, long ind, long prec)
    1019             : {
    1020             :   GEN res;
    1021       24486 :   switch (labs(ind))
    1022             :   {
    1023       12572 :     case 1: res = F21taylor1(a, b, c, z, prec); break;
    1024       11242 :     case 2: res = F21taylor(a, gsub(c, b), c, gdiv(z, gsubgs(z, 1)), prec);
    1025       11242 :             res = gmul(res, gpow(gsubsg(1,z), gneg(a), prec)); break;
    1026         336 :     case 3: res = F21taylor(a, b, c, z, prec); break;
    1027         168 :     case 4: res = F21taylor4(a, b, c, z, prec); break;
    1028          84 :     case 5: res = F21taylor5(a, b, c, z, prec); break;
    1029          84 :     default:res = F21taylor6(a, b, c, z, prec); break;
    1030             :   }
    1031       24486 :   return gprec_wtrunc(res, prec);
    1032             : }
    1033             : 
    1034             : static long
    1035        3108 : hypersimplify(GEN *pn, GEN *pd)
    1036             : {
    1037        3108 :   GEN n = *pn, d = *pd;
    1038        3108 :   long j, ld = lg(d), ln = lg(n);
    1039        5922 :   for (j = 1; j < ld; j++)
    1040             :   {
    1041        3066 :     GEN t = gel(d, j);
    1042             :     long k;
    1043        7777 :     for (k = 1; k < ln; k++)
    1044        4963 :       if (gequal(t, gel(n, k)))
    1045             :       {
    1046         252 :         *pd = vecsplice(d, j);
    1047         252 :         *pn = vecsplice(n, k); return hypersimplify(pd, pn) + 1;
    1048             :       }
    1049             :   }
    1050        2856 :   return 0;
    1051             : }
    1052             : 
    1053             : static GEN
    1054        2086 : f_pochall(void *E, GEN n)
    1055             : {
    1056        2086 :   GEN S, a, tmp, N = gel(E,1), D = gel(E,2);
    1057        2086 :   long j, prec = itou(gel(E,3));
    1058        2086 :   S = gen_0;
    1059        9471 :   for (j = 1; j < lg(N); j++)
    1060             :   {
    1061        7385 :     a = gel(N, j);
    1062        7385 :     tmp = gsub(glngamma(gadd(n, a), prec), glngamma(a, prec));
    1063        7385 :     S = gadd(S, tmp);
    1064             :   }
    1065        7385 :   for (j = 1; j < lg(D); j++)
    1066             :   {
    1067        5299 :     a = gel(D, j);
    1068        5299 :     tmp = gsub(glngamma(gadd(n, a), prec), glngamma(a, prec));
    1069        5299 :     S = gsub(S, tmp);
    1070             :   }
    1071        2086 :   return gexp(gsub(S, glngamma(gaddsg(1, n), prec)), prec);
    1072             : }
    1073             : static GEN
    1074         371 : f_pochall_alt(void *E, GEN n)
    1075         371 : { GEN z = f_pochall(E,n); return mpodd(n)? gneg(z): z; }
    1076             : /* z = \pm1 */
    1077             : static GEN
    1078          63 : sumz(GEN N, GEN D, long z, long prec)
    1079             : {
    1080          63 :   void *E = (void*)mkvec3(N, D, utoi(prec));
    1081             :   GEN tab, be;
    1082          63 :   if (z == -1) return sumalt(E, f_pochall_alt, gen_0, prec);
    1083          56 :   be = gsub(vecsum(D), vecsum(N)); check_hyp1(be);
    1084          56 :   tab = sumnummonieninit(be, NULL, gen_0, prec);
    1085          56 :   return sumnummonien(E, f_pochall, gen_0, tab, prec);
    1086             : }
    1087             : 
    1088             : static GEN
    1089        5712 : hypergeom_arg(GEN x)
    1090             : {
    1091        5712 :   if (!x) return cgetg(1,t_VEC);
    1092        5712 :   return (typ(x) == t_VEC)? x: mkvec(x);
    1093             : }
    1094             : 
    1095             : static GEN
    1096        3150 : hypergeom_i(GEN N, GEN D, GEN z, long prec)
    1097             : {
    1098             :   long nN, nD;
    1099        3150 :   if (!is_scalar_t(typ(z))) pari_err_TYPE("hypergeom",z);
    1100        3150 :   if (gequal0(z)) return gen_1;
    1101        3150 :   nN = lg(N) - 1;
    1102        3150 :   nD = lg(D) - 1;
    1103        3150 :   if (nD >= (nN? nN: 2)) return Ftaylor(N, D, z, prec);
    1104        2100 :   if (nD == nN - 1 && nN >= 3)
    1105             :   {
    1106         126 :     GEN d = gsubsg(1, gabs(z,LOWDEFAULTPREC));
    1107         126 :     long ed = gexpo(d);
    1108             :     /* z in unit disc but "away" from unit circle */
    1109         126 :     if (gsigne(d) > 0 && ed > -prec2nbits(prec)/4
    1110          21 :         && (nN != 3 || ed > -15)) /* For 3F2 we can use integral */
    1111          14 :       return Ftaylor(N, D, z, prec);
    1112         112 :     if (gequal1(z))  return sumz(N, D, 1, prec);
    1113          56 :     if (gequalm1(z)) return sumz(N, D,-1, prec);
    1114             :   }
    1115        2023 :   switch (nN)
    1116             :   {
    1117         112 :     case 0:
    1118         112 :       if (nD == 0) return gexp(z, prec);
    1119          21 :       if (nD == 1) return F01(gel(D,1), z, prec);
    1120         350 :     case 1: return F10(gel(N, 1), z, prec);
    1121        1449 :     case 2:
    1122        1449 :       if (nD == 0) return F20(gel(N,1), gel(N,2), z, prec);
    1123        1393 :       if (nD == 1) return F21(gel(N,1), gel(N,2), gel(D,1), z, prec);
    1124             :     case 3:
    1125         105 :       if (nD == 0) break;
    1126         105 :       if (nD == 1) return F31(gel(N,1), gel(N,2), gel(N,3), gel(D,1), z, prec);
    1127          49 :       if (nD == 2) return F32(N, D, z, prec);
    1128             :   }
    1129           7 :   pari_err_IMPL("this hypergeometric function");
    1130             :   return NULL; /*LCOV_EXCL_LINE*/
    1131             : }
    1132             : 
    1133             : /* [v[i]+k, i=1..#v]; no need to optimize for k = 0 */
    1134             : static GEN
    1135         728 : RgV_z_add(GEN v, long k)
    1136             : {
    1137             :   long i, l;
    1138         728 :   GEN w = cgetg_copy(v, &l);
    1139        1456 :   for (i = 1; i < l; i++) gel(w,i) = gaddgs(gel(v,i), k);
    1140         728 :   return w;
    1141             : }
    1142             : static GEN
    1143        1442 : vpoch_mul(GEN v, long k)
    1144             : {
    1145        1442 :   long i, l = lg(v);
    1146             :   GEN p;
    1147        1442 :   if (l == 1) return gen_1;
    1148         903 :   p = gmael(v, 1, k);
    1149        1204 :   for (i = 2; i < l; i++) p = gmul(p, gmael(v, i, k));
    1150         903 :   return p;
    1151             : }
    1152             : 
    1153             : static GEN
    1154          56 : serhypergeom(GEN N, GEN D, GEN y, long prec)
    1155             : {
    1156             :   pari_sp av;
    1157          56 :   GEN S, pn, vN, vD, y0 = NULL;
    1158             :   long v, l, i;
    1159          56 :   if (!signe(y)) return gadd(gen_1, y);
    1160          49 :   v = valp(y);
    1161          49 :   if (v < 0) pari_err_DOMAIN("hypergeom","valuation", "<", gen_0, y);
    1162          49 :   l = lg(y);
    1163          49 :   if (v) S = gen_1;
    1164             :   else
    1165             :   {
    1166          28 :     y0 = gel(y, 2);
    1167          28 :     y = serchop0(y);
    1168          28 :     l = 3 + (l - 3) / valp(y);
    1169          28 :     S = hypergeom(N, D, y0, prec);
    1170             :   }
    1171          49 :   vN = RgV_vpoch(N, l-1);
    1172          49 :   vD = RgV_vpoch(D, l-1); av = avma;
    1173          49 :   pn = y;
    1174         770 :   for (i = 1; i < l; i++)
    1175             :   {
    1176         721 :     GEN H = gdiv(vpoch_mul(vN, i), vpoch_mul(vD, i));
    1177         721 :     if (y0) H = gmul(H, hypergeom_i(RgV_z_add(N,i), RgV_z_add(D,i), y0,prec));
    1178         721 :     S = gadd(S, gmul(pn, H)); if (i + 1 < l) pn = gdivgu(gmul(pn, y), i + 1);
    1179         721 :     if (gc_needed(av,1))
    1180             :     {
    1181           0 :       if (DEBUGMEM>1) pari_warn(warnmem,"hypergeom, i = %ld / %ld", i,l-1);
    1182           0 :       gerepileall(av, 2, &S, &pn);
    1183             :     }
    1184             :   }
    1185          49 :   return S;
    1186             : }
    1187             : 
    1188             : GEN
    1189        2856 : hypergeom(GEN N, GEN D, GEN y, long prec)
    1190             : {
    1191        2856 :   pari_sp av = avma;
    1192             :   GEN z;
    1193             :   long j, n;
    1194        2856 :   N = hypergeom_arg(N);
    1195        2856 :   D = hypergeom_arg(D);
    1196        2856 :   n = hypersimplify(&N, &D);
    1197        5516 :   for (j = 1; j < lg(D); j++)
    1198        2674 :     if (isnegint(gel(D,j)))
    1199          14 :       pari_err_DOMAIN("hypergeom", stack_sprintf("b[%ld]", j + n),
    1200          14 :                       "<=", gen_0, gel(D,j));
    1201        2842 :   if (is_scalar_t(typ(y)))
    1202        2786 :     return gerepilecopy(av, hypergeom_i(N, D, y, prec));
    1203          56 :   if (!(z = toser_i(y))) pari_err_TYPE("hypergeom", y);
    1204          56 :   return gerepileupto(av, serhypergeom(N, D, z, prec));
    1205             : }

Generated by: LCOV version 1.13