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

Generated by: LCOV version 1.13