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

Generated by: LCOV version 1.14