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 116039 : isnegint2(GEN a, long *pa)
30 : {
31 : GEN b;
32 116039 : if (!is_scalar_t(typ(a)) || !gequal0(imag_i(a))) return 0;
33 94122 : a = real_i(a); if (gsigne(a) > 0) return 0;
34 6076 : b = ground(a); if (!gequal(a, b)) return 0;
35 4991 : if (pa) *pa = -itos(b);
36 4991 : return 1;
37 : }
38 : static int
39 27825 : isnegint(GEN a) { return isnegint2(a, NULL); }
40 : static int
41 58884 : isnegint_approx(GEN a, long bit)
42 : {
43 : GEN b;
44 58884 : if (gexpo(imag_i(a)) > -bit) return 0;
45 48076 : a = real_i(a); if (gsigne(a) > 0) return 0;
46 18375 : b = ground(a); return gexpo(gsub(a, b)) < -bit;
47 : }
48 : static int
49 73885 : is0(GEN a, long bit) { return gequal0(a) || gexpo(a) < -bit; }
50 : static int
51 12936 : islong(GEN a, long *m, long prec)
52 : {
53 12936 : *m = itos(ground(real_i(a)));
54 12936 : 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 gc_GEN(av, airy_i(z, prec)); }
105 :
106 : /* Gamma(a)*Gamma(b); allow a = NULL [omit] */
107 : static GEN
108 28819 : mulgamma2(GEN a, GEN b, long prec)
109 28819 : { GEN gb = ggamma(b, prec); return a? gmul(ggamma(a, prec), gb): gb; }
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 28602 : mulgammav2(GEN v, long prec)
117 28602 : { 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 = precdbl(l);
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)) (void)gc_all(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)) (void)gc_all(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 : (void)gc_all(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 gc_GEN(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 gc_upto(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 gc_upto(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 114352 : poch(GEN a, long n, long prec)
345 : {
346 114352 : GEN S = real_1(prec);
347 : long j;
348 1387680 : for (j = 0; j < n; j++) S = gmul(S, gaddsg(j, a));
349 114352 : 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 114086 : Npoch(GEN a, long n) { return gnorm(poch(a, n, LOWDEFAULTPREC)); }
363 : static GEN
364 38472 : Npochden(GEN a, long n)
365 : {
366 38472 : GEN L = Npoch(a, n), r = ground(real_i(a));
367 38472 : 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 38472 : return L;
373 : }
374 :
375 : /* |x + z|^2 */
376 : static GEN
377 111909 : normpol2(GEN z)
378 : {
379 111909 : GEN a = real_i(z), b = imag_i(z);
380 111909 : 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 75362 : vnormpol2(GEN v)
385 : {
386 75362 : long i, l = lg(v);
387 : GEN P;
388 75362 : if (l == 1) return pol_1(0);
389 75362 : P = normpol2(gel(v,1));
390 111909 : for (i = 2; i < l; i++) P = RgX_mul(P, normpol2(gel(v,i)));
391 75362 : return P;
392 : }
393 :
394 : /* return an 'extraprec' */
395 : static long
396 37681 : precFtaylor(GEN N, GEN D, GEN z, long *pmi)
397 : {
398 37681 : GEN v, ma, P = vnormpol2(D), Q = vnormpol2(N), Nz = gnorm(z);
399 37681 : double wma, logNz = (gexpo(Nz) < -27)? -27: dbllog2(Nz) / 2;
400 37681 : long pr = LOWDEFAULTPREC, prec = precision(z);
401 37681 : long lN = lg(N), lD = lg(D), j, i, lv;
402 :
403 37681 : P = RgX_shift_shallow(P, 2);
404 : /* avoid almost cancellation of leading coeff if |z| ~ 1 */
405 37681 : if (!prec || fabs(logNz) > 1e-38) Q = RgX_Rg_mul(Q,Nz);
406 111881 : for (j = 1, ma = NULL; j < lN; j++)
407 : {
408 74200 : GEN Nj = gel(N,j);
409 74200 : if (isint(Nj,&Nj) && signe(Nj) <= 0
410 5236 : && (!ma || abscmpii(ma, Nj) < 0)) ma = Nj;
411 : }
412 : /* use less sharp fujiwara_bound_real(,1) ? */
413 37681 : v = realroots(gsub(P,Q), mkvec2(gen_0,mkoo()), pr);
414 37681 : v = ZV_to_zv(ground(v)); lv = lg(v);
415 37681 : if (ma)
416 : {
417 5061 : long sma = is_bigint(ma)? LONG_MAX: maxss(labs(itos(ma))-1, 1);
418 10773 : for (i = 1; i < lv; i++) v[i] = maxss(minss(sma, v[i]), 1);
419 : }
420 76125 : for (i = 1, wma = 0.; i < lv; i++)
421 : {
422 38444 : GEN t1 = gen_1, t2 = gen_1;
423 38444 : long u = v[i];
424 : double t;
425 114058 : for (j = 1; j < lN; j++) t1 = gmul(t1, Npoch(gel(N,j), u));
426 76916 : for (j = 1; j < lD; j++) t2 = gmul(t2, Npochden(gel(D,j), u));
427 38444 : t = dbllog2(gdiv(t1,t2)) / 2 + u * logNz - dbllog2(mpfactr(u,pr));
428 38444 : wma = maxdd(wma, t); /* t ~ log2 | N(u)/D(u) z^u/u! | */
429 : }
430 : /* make up for exponential decrease in exp() */
431 37681 : if (gsigne(real_i(z)) < 0) wma -= gtodouble(real_i(z)) / M_LN2;
432 37681 : *pmi = maxss(v[lv-1],1); return nbits2extraprec(wma+BITS_IN_LONG);
433 : }
434 :
435 : static GEN
436 24430 : Ftaylor(GEN N, GEN D, GEN z, long prec)
437 : {
438 : pari_sp av;
439 : GEN C, S;
440 24430 : long i, j, ct, lN = lg(N), lD = lg(D), pradd, mi, bitmin, tol;
441 24430 : pradd = precFtaylor(N, D, z, &mi);
442 24430 : if (pradd > 0)
443 : {
444 24430 : prec += pradd;
445 24430 : N = gprec_wensure(N, prec);
446 24430 : D = gprec_wensure(D, prec);
447 24430 : z = gprec_wensure(z, prec);
448 : }
449 24430 : bitmin = -(prec2nbits(prec) + 10);
450 24430 : S = C = real_1(prec); ct = 0; j = 0; tol = 0;
451 24430 : av = avma;
452 : for(;;)
453 7319436 : {
454 7343866 : GEN a = gen_1, b = gen_1;
455 21360848 : for (i = 1; i < lN; i++) a = gmul(a, gaddsg(j, gel(N,i)));
456 14688715 : for (i = 1; i < lD; i++) b = gmul(b, gaddsg(j, gel(D,i)));
457 7343866 : C = gmul(C, gmul(gdiv(a, gmulsg(j+1, b)), z));
458 7343866 : if (gequal0(C)) break;
459 7343586 : if (j > mi) tol = gequal0(S)? 0: gexpo(C) - gexpo(S);
460 7343586 : S = gadd(S, C); j++;
461 7343586 : if (j > mi)
462 7114063 : { if (tol > bitmin) ct = 0; else if (++ct >= lN+lD-2) break; }
463 7319436 : if (gc_needed(av, 1)) (void)gc_all(av, 2, &S, &C);
464 : }
465 24430 : return S;
466 : }
467 :
468 : static GEN
469 210 : bind(GEN a, GEN b, GEN c, long ind)
470 : {
471 210 : switch(ind)
472 : {
473 119 : case 1: case 2: return gsub(c, b);
474 63 : case 5: case 6: return gsub(gaddsg(1, a), c);
475 28 : default: return b; /* 3,4 */
476 : }
477 : }
478 : static GEN
479 59094 : cind(GEN a, GEN b, GEN c, long ind)
480 : {
481 59094 : switch(ind)
482 : {
483 29505 : case 1: case 6: return gaddsg(1, gsub(a,b));
484 29498 : case 4: case 5: return gsub(gaddsg(1, gadd(a,b)), c);
485 91 : default: return c; /* 2,3 */
486 : }
487 : }
488 : static GEN
489 111524 : zind(GEN z, long ind)
490 : {
491 111524 : switch(ind)
492 : {
493 24850 : case 1: return ginv(gsubsg(1, z));
494 29533 : case 2: return gdiv(z, gsubgs(z, 1));
495 16142 : case 4: return gsubsg(1, z);
496 16142 : case 5: return gsubsg(1, ginv(z));
497 24857 : case 6: return ginv(z);
498 : default: return z; /* 3, LCOV_EXCL_LINE*/
499 : }
500 : }
501 :
502 : /* z not 0 or 1, c not a nonpositive integer. Among the generic 6
503 : * Kummer transforms 1/(1-z),z/(z-1),z,1-z,1-1/z,1/z, check in this order
504 : * which ones are allowed and return the index for the one with smallest
505 : * resulting |z|. If that modulus is "far enough" from 1 [we can sum the Taylor
506 : * series], return it negated. Otherwise, we'll use integral representations. */
507 : static long
508 29442 : F21ind(GEN a, GEN b, GEN c, GEN z, long bit)
509 : {
510 29442 : GEN v = const_vec(6, mkoo());
511 29442 : long ind = 0, B = bit - 5;
512 29442 : const long LD = LOWDEFAULTPREC;
513 29442 : if (!isnegint_approx(cind(a,b,c, 1),B))
514 : { /* 1 and 6 are allowed */
515 24822 : gel(v,1) = gabs(zind(z,1), LD);
516 24822 : gel(v,6) = gabs(zind(z,6), LD);
517 : }
518 29442 : gel(v,2) = gabs(zind(z,2), LD);
519 29442 : gel(v,3) = gabs(z, LD);
520 29442 : if (!isnegint_approx(cind(a,b,c, 4),B))
521 : { /* 4 and 5 are allowed */
522 16114 : gel(v,4) = gabs(zind(z,4), LD);
523 16114 : gel(v,5) = gabs(zind(z,5), LD);
524 : }
525 29442 : ind = vecindexmin(v); /* |znew| <= 1; close to 1 ? */
526 29442 : return (gexpo(gsubgs(gel(v,ind),1)) > -maxss(bit / 4, 32))? -ind: ind;
527 : }
528 : static GEN
529 73738 : mul3(GEN a, GEN b, GEN c) { return gmul(a,gmul(b, c)); }
530 :
531 : /* return the polynomial of degree N >= 1
532 : * \sum_{i=0} binomial(N,i) (c-a)_i (a)_{N-i} X^i */
533 : static GEN
534 28 : F21lam(long N, GEN a, GEN c)
535 : {
536 28 : GEN C = vecbinomial(N), S = cgetg(N+2, t_VEC);
537 28 : GEN vb = vpoch(gsub(c,a), N), va = vpoch(a, N);
538 : long i;
539 28 : gel(S,1) = gel(va,N);
540 28 : for (i = 1; i < N; i++) gel(S,i+1) = mul3(gel(C,i+1), gel(vb,i), gel(va,N-i));
541 28 : gel(S,i+1) = gel(vb,N); return RgV_to_RgX(S,0);
542 : }
543 :
544 : /* F(-m,b; c; z), m >= 0; c != 0, -1, -2, ... */
545 : static GEN
546 4746 : F21finitetaylor(long m, GEN b, GEN c, GEN z, long prec)
547 : {
548 : pari_sp av;
549 : GEN C, P, S;
550 : long j, ct, pradd, mi, bitmin;
551 4746 : if (isnegint2(b, &j) && j < m) { b = stoi(-m); m = j; }
552 4746 : pradd = precFtaylor(mkvec2(stoi(-m), b), mkvec(c), z, &mi);
553 4746 : if (pradd > 0)
554 : {
555 4746 : prec += pradd;
556 4746 : b = gprec_wensure(b, prec);
557 4746 : c = gprec_wensure(c, prec);
558 4746 : z = gprec_wensure(z, prec);
559 : }
560 4746 : bitmin = -(prec2nbits(prec) + 10);
561 4746 : C = vecbinomial(m);
562 4746 : P = real_1(prec); S = P; ct = 0;
563 4746 : av = avma;
564 133189 : for(j = 1; j <= m; j++) /* NIST 15.2.4 */
565 : {
566 : GEN p;
567 128443 : P = gmul(P, gmul(z, gdiv(gaddsg(j-1,b), gaddsg(j-1,c)))); /*(b)_j z^j / (c)_j*/
568 128443 : p = gmul(P, gel(C, j+1)); /* times binomial(m,j) */
569 128443 : if (j > mi && !gequal0(S))
570 7336 : { if (gexpo(p) - gexpo(S) > bitmin) ct = 0; else if (++ct == 3) break; }
571 128443 : S = odd(j)? gsub(S, p): gadd(S, p);
572 128443 : if (gc_needed(av, 1)) (void)gc_all(av, 2, &S, &P);
573 : }
574 4746 : return S;
575 : }
576 :
577 : static GEN
578 119 : F21finite_i(GEN F, GEN B, GEN C, long m, long prec)
579 119 : { return gdiv(gmul(F, poch(B, m, prec)), poch(C, m, prec)); }
580 :
581 : /* F(-m,b; c; z), m >= 0; c not a nonpositive integer */
582 : static GEN
583 4746 : F21finite(long m, GEN b, GEN c, GEN z, long prec)
584 : {
585 4746 : GEN a = stoi(-m), b1, c1, z1, F;
586 4746 : long ind = F21ind(a, b, c, z, prec2nbits(prec));
587 :
588 4746 : if (ind > 0 || ind == -3) return F21finitetaylor(m, b, c, z, prec);
589 210 : ind = labs(ind);
590 210 : z1 = zind(z, ind);
591 210 : b1 = bind(a, b, c, ind);
592 210 : c1 = cind(a, b, c, ind); /* not a nonpositive integer */
593 210 : F = F21finitetaylor(m, b1, c1, z1, prec);
594 210 : switch (ind)
595 : {
596 28 : case 1: return F21finite_i(gmul(F, gpowgs(gsubsg(1,z), m)), b, c, m, prec);
597 91 : case 2: return gmul(F, gpowgs(gsubsg(1,z), m));
598 28 : case 4: return F21finite_i(F, gsub(c,b), c, m, prec);
599 28 : case 5: return F21finite_i(gmul(F, gpowgs(z,m)), gsub(c,b), c, m, prec);
600 35 : case 6: return F21finite_i(gmul(F, gpowgs(gneg(z),m)), b, c, m, prec);
601 : }
602 : return NULL; /*LCOV_EXCL_LINE*/
603 : }
604 :
605 : /**********************************************************/
606 :
607 : /* allow a = NULL [omit] */
608 : static GEN
609 28 : multgam(GEN a, GEN b, GEN c, GEN d, long prec)
610 : {
611 28 : if (isnegint(c) || isnegint(d)) return gen_0;
612 28 : return gdiv(mulgamma2(a, b, prec), mulgamma2(c, d, prec));
613 : }
614 :
615 : static GEN
616 119 : intnumsplit(void *E, GEN (*f)(void*, GEN), GEN a, GEN b, GEN z, long prec)
617 : {
618 119 : if (!z) return intnum(E, f, a, b, NULL, prec);
619 7 : return gadd(intnum(E, f, a, z, NULL, prec),
620 : intnum(E, f, z, b, NULL, prec));
621 : }
622 :
623 : /* (1 - zt)^a T(1-zt) t^b (1-t)^c */
624 : static GEN
625 73738 : fF21(void *E, GEN t)
626 : {
627 73738 : GEN z = gel(E,1), a = gel(E,2), b = gel(E,3), c = gel(E,4), T = gel(E,5);
628 73738 : GEN u = gsubsg(1, gmul(z, t));
629 73738 : long prec = precision(t);
630 73738 : GEN x = mul3(gpow(u, a, prec), gpow(t, b, prec), gpow(gsubsg(1,t), c, prec));
631 73738 : return T? gmul(x, poleval(T, u)): x;
632 : }
633 :
634 : /* z != 1 */
635 : static GEN
636 119 : int21(GEN a, GEN b, GEN c, GEN z, GEN T, long prec)
637 : {
638 119 : GEN A = gneg(a), B = gsubgs(b,1), C = gsubgs(gsub(c,b), 1), rA = real_i(A);
639 119 : GEN pz = NULL, p0 = mkendpt(gen_0, B), p1 = mkendpt(gen_1, C);
640 119 : GEN E = mkvec5(z, A, B, C, T);
641 :
642 119 : if (gcmpgs(rA, 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,rA);
649 112 : else if (gcmpgs(rA,-1) <= 0) prec += ((gexpo(rA)+1)>>1) * EXTRAPREC64;
650 119 : return intnumsplit((void*)E, fF21, p0, p1, pz, prec);
651 : }
652 :
653 : static GEN F21taylorind(GEN a, GEN b, GEN c, GEN z, long ind, long prec);
654 : /* c not a nonpositive integer */
655 : static GEN
656 30331 : F21_i(GEN a, GEN b, GEN c, GEN z, long prec)
657 : {
658 30331 : GEN res, p = NULL, T = NULL;
659 30331 : long m, ind, prec2, bitprec = prec2nbits(prec);
660 30331 : if (is0(imag_i(z), bitprec)) z = real_i(z);
661 30331 : if (is0(z, bitprec)) return real_1(prec);
662 29477 : if (gequal1(z))
663 : {
664 35 : GEN x = gsub(c, gadd(a, b)); check_hyp1(x);
665 28 : return multgam(c, x, gsub(c,a), gsub(c,b), prec);
666 : }
667 : /* z != 0,1 */
668 29442 : if (isnegint2(b, &m)) return F21finite(m, a, c, z, prec);
669 29239 : if (isnegint2(a, &m)) return F21finite(m, b, c, z, prec);
670 24787 : if (isnegint(gsub(c, b))) swap(a, b);
671 24787 : if (isnegint2(gsub(c, a), &m))
672 : {
673 91 : GEN x = gpow(gsubsg(1, z), gneg(gaddsg(m, b)), prec);
674 91 : return gmul(x, F21finite(m, gsub(c, b), c, z, prec));
675 : }
676 : /* None of a, b, c, c-a, c-b is a nonpositive integer.
677 : * Try Kummer transforms */
678 24696 : ind = F21ind(a, b, c, z, bitprec);
679 24696 : prec2 = prec + EXTRAPREC64;
680 24696 : a = gprec_wensure(a,prec2);
681 24696 : b = gprec_wensure(b,prec2);
682 24696 : c = gprec_wensure(c,prec2);
683 24696 : z = gprec_wensure(z,prec2);
684 : /* can ensure |z| << 1: sum the Taylor series */
685 24696 : if (ind < 0) return gprec_wtrunc(F21taylorind(a,b,c, z, ind, prec), prec);
686 : /* can ensure |z| <= 1, but too close to 1: use integral 15.6.1 for
687 : * F21 / Gamma(c) */
688 126 : if (gsigne(real_i(b)) <= 0)
689 : {
690 35 : if (gsigne(real_i(a)) <= 0)
691 : { /* 15.5.16 (a and b swapped) divided by c */
692 : GEN p1,p2;
693 7 : if (gcmp(real_i(b), real_i(a)) < 0) swap(a,b);
694 : /* FIXME: solve recursion as below with F21lam */
695 7 : p1 = gmul(gsubsg(1, z), F21_i(a, gaddsg(1,b), c, z, prec));
696 7 : p2 = gmul(gmul(gsubsg(1, gdiv(a,c)), z),
697 : F21_i(a, gaddsg(1,b), gaddsg(1,c), z, prec));
698 7 : return gadd(p1, p2);
699 : }
700 28 : swap(a,b);
701 : }
702 : /* Here real(b) > 0 */
703 119 : if (gcmp(real_i(a), real_i(b)) < 0 && gsigne(real_i(a)) > 0) swap(a,b);
704 : /* further, either real(a) <= 0 or real(a) >= real(b) */
705 119 : if (gcmp(real_i(c), real_i(b)) <= 0)
706 : { /* T(x) = sum_i binomial(N,i) (c-a)_i (a)_{N-i} x^i
707 : * (c)_N F(a,b,c;z) = sum_{i = 0}^N Ti F(a+N-i,b,c+N;z) for N >= 0,
708 : * where F = 2F1 / Gamma(c) */
709 28 : long N = 1 + itos(gfloor(gsub(real_i(b),real_i(c)))); /* >= 1 */
710 28 : T = F21lam(N, a, c);
711 28 : p = poch(c, N, prec);
712 28 : c = gaddsg(N,c); /* now real(b) < real(c) */
713 28 : a = gaddsg(N,a);
714 : }
715 : /* real(b) < real(c) */
716 119 : res = int21(a, b, c, z, T, prec2); if (p) res = gdiv(res, p);
717 : /* res = 2F1(a0,b0,c0; z)*Gamma(b)Gamma(c-b)/Gamma(c) by 15.6.1 */
718 119 : return gdiv(gmul(ggamma(c,prec),res), mulgamma2(b, gsub(c,b), prec));
719 : }
720 :
721 : /* c not a non-positive integer */
722 : static GEN
723 30261 : F21(GEN a, GEN b, GEN c, GEN z, long prec)
724 : {
725 30261 : GEN res = F21_i(a, b, c, z, prec);
726 30254 : long ex = labs(gexpo(res)), bitprec = prec2nbits(prec);
727 30254 : if (ex > bitprec)
728 : {
729 56 : prec = nbits2prec(ex + bitprec);
730 56 : res = F21_i(gprec_wensure(a,prec), gprec_wensure(b,prec),
731 : gprec_wensure(c,prec), gprec_wensure(z,prec), prec);
732 : }
733 30254 : return res;
734 : }
735 :
736 : static GEN
737 11592 : F21taylor(GEN a, GEN b, GEN c, GEN z, long prec)
738 : {
739 11592 : pari_sp av = avma;
740 11592 : GEN r = gdiv(Ftaylor(mkvec2(a,b), mkvec(c), z, prec), ggamma(c, prec));
741 11592 : return gc_upto(av, r);
742 : }
743 :
744 : static GEN
745 8505 : F21taylorlim(GEN N, long m, GEN z, GEN Z, long ind, long prec)
746 : {
747 : pari_sp av;
748 : GEN C, P, S, a, b;
749 8505 : long j, ct, pradd, mi, fl, bitmin, tol, si = (ind == 5 || ind == 6)? -1: 1;
750 8505 : pradd = precFtaylor(N, mkvec(stoi(m + 1)), z, &mi);
751 8505 : if (pradd)
752 : {
753 8505 : prec += pradd;
754 8505 : N = gprec_wensure(N, prec);
755 8505 : z = gprec_wensure(z, prec);
756 8505 : Z = gprec_wensure(Z, prec);
757 : }
758 8505 : av = avma; a = gel(N,1); b = gel(N,2);
759 8505 : bitmin = -(prec2nbits(prec) + 10);
760 8505 : P = glog(Z, prec); if (ind == 4 || ind == 5) P = gneg(P);
761 8505 : P = gadd(P, gsub(gpsi(stoi(m+1), prec), mpeuler(prec)));
762 8505 : P = gsub(P, gadd(gpsi(a, prec), gpsi(si == -1? gsubsg(1, b): b, prec)));
763 8505 : C = real_1(prec); ct = 0; tol = 0; S = P;
764 8505 : for(j = 0, fl = 1;;)
765 994800 : {
766 1003305 : GEN v1 = gaddgs(a, j), v2 = gaddgs(b, j);
767 1003305 : long J = (j+1) * (j+1+m);
768 1003305 : C = gmul(C, gdivgs(gmul(z, v1), J));
769 1003305 : if (gequal0(v2)) fl = 0; else C = gmul(C, v2);
770 1003305 : if (j > mi) tol = gequal0(S) ? 0 : gexpo(C) - gexpo(S);
771 1003305 : if (fl)
772 : {
773 1002027 : P = gadd(P, gsub(uutoQ(2*j+2+m, J), gadd(ginv(v1), ginv(v2))));
774 1002027 : S = gadd(S, gmul(C, P));
775 : }
776 : else
777 1278 : S = (si == 1)? gadd(S, C): gsub(S, C);
778 1003305 : if (++j > mi) { if (tol > bitmin) ct = 0; else if (++ct == 3) break; }
779 994800 : if (gc_needed(av, 1)) (void)gc_all(av, 3, &S, &C, &P);
780 : }
781 8505 : return gdiv(S, mpfact(m));
782 : }
783 :
784 : /* N = [a,b]; (m-1)! sum_{0 <= k < m} (-z)^k (a)_k (b)_k / k!(m-k)...(m-1) */
785 : static GEN
786 8505 : F21finitelim(GEN N, long m, GEN z, long prec)
787 : {
788 : GEN C, S, a, b;
789 : long k;
790 8505 : if (!m) return gen_0;
791 112 : a = gel(N,1);
792 112 : b = gel(N,2);
793 112 : S = C = real_1(prec + EXTRAPREC64);
794 175 : for (k = 1; k < m; k++)
795 : {
796 63 : GEN v1 = gaddsg(k-1, a), v2 = gaddsg(k-1, b);
797 63 : C = gdivgs(gmul(C, gmul(gmul(v1, v2), z)), k*(k-m));
798 63 : S = gadd(S, C);
799 : }
800 112 : return gmul(S, mpfact(m-1));
801 : }
802 :
803 : static GEN
804 5796 : OK_gadd(GEN x, GEN y, long prec0, long *pprec,
805 : GEN *d1,GEN *d2,GEN *d3,GEN *d4,GEN *d5,GEN *d6,GEN *d7,GEN *d8)
806 : {
807 5796 : long prec = *pprec;
808 5796 : GEN z = gadd(x,y);
809 5796 : if (!gequal0(z) && gexpo(z)-gexpo(x) >= prec2nbits(prec0) - prec2nbits(prec))
810 4431 : return z;
811 1365 : *pprec = prec = precdbl(prec);
812 1365 : *d1 = gprec_wensure(*d1, prec); *d2 = gprec_wensure(*d2, prec);
813 1365 : *d3 = gprec_wensure(*d3, prec); *d4 = gprec_wensure(*d4, prec);
814 1365 : *d5 = gprec_wensure(*d5, prec); *d6 = gprec_wensure(*d6, prec);
815 1365 : *d7 = gprec_wensure(*d7, prec); *d8 = gprec_wensure(*d8, prec);
816 1365 : return NULL;
817 : }
818 :
819 : static GEN
820 4431 : FBaux1(GEN v1, GEN g1, GEN c1, GEN v2, GEN g2, GEN c2, GEN z, GEN bma,
821 : long prec0, long prec)
822 : {
823 4431 : GEN pi = mppi(prec);
824 : for (;;)
825 1365 : {
826 5796 : GEN t1 = gdiv(c1, mulgammav2(g1, prec)), r1;
827 5796 : GEN t2 = gdiv(c2, mulgammav2(g2, prec)), r2, F;
828 5796 : r1 = gmul(t1, F21taylor(gel(v1,1), gel(v1,2), gel(v1,3), z, prec));
829 5796 : r2 = gmul(t2, F21taylor(gel(v2,1), gel(v2,2), gel(v2,3), z, prec));
830 5796 : F = OK_gadd(r1, r2, prec0, &prec, &c1,&c2, &g1,&g2, &v1,&v2, &z,&bma);
831 5796 : if (F) return gmul(F, gdiv(pi, gsin(gmul(pi, bma), prec)));
832 : }
833 : }
834 :
835 : static GEN
836 8505 : FBaux2(GEN v1, GEN g1, GEN c1, long m, GEN z1, GEN c2, GEN g2, GEN v2, GEN z2,
837 : GEN Z2, long ind, long prec)
838 : {
839 8505 : GEN t1 = gdiv(c1, mulgammav2(g1, prec)), r1;
840 8505 : GEN t2 = gdiv(c2, mulgammav2(g2, prec)), r2;
841 8505 : r1 = gmul(t1, F21finitelim(v1, m, z1, prec));
842 8505 : r2 = gmul(t2, F21taylorlim(v2, m, z2, Z2, ind, prec));
843 8505 : return gadd(r1, r2);
844 : }
845 :
846 : /* 1 / (1-z) */
847 : static GEN
848 12572 : F21taylor1(GEN a, GEN b, GEN c, GEN z, long prec)
849 : {
850 12572 : GEN bma = gsub(b, a), coe1, coe2, z1, g1, g2, v1, v2, Z;
851 : long m;
852 12572 : if (!islong(bma,&m, prec))
853 : {
854 : GEN b1, c1, e1, c2;
855 4207 : b1 = gsub(c, b);
856 4207 : c1 = gsubsg(1, bma);
857 4207 : e1 = gsub(c, a);
858 4207 : coe1 = gpow(gsubsg(1, z), gneg(a), prec);
859 4207 : c2 = gaddgs(bma, 1);
860 4207 : coe2 = gneg(gpow(gsubsg(1, z), gneg(b), prec));
861 4207 : z1 = ginv(gsubsg(1, z));
862 4207 : return FBaux1(mkvec3(a,b1,c1), mkvec2(b,e1), coe1, mkvec3(b,e1,c2),
863 : mkvec2(a,b1), coe2, z1, bma, prec, prec);
864 : }
865 8365 : if (m < 0) { swap(a,b); m = -m; }
866 8365 : Z = gsubsg(1, z);
867 8365 : coe1 = gpow(Z, gneg(a), prec);
868 8365 : v2 = g1 = mkvec2(gaddgs(a, m), gsub(c, a));
869 8365 : v1 = mkvec2(a, gsub(c, gaddgs(a, m)));
870 8365 : z1 = ginv(Z);
871 8365 : coe2 = gmul(coe1, gpowgs(z1, m)); if (m & 1L) coe2 = gneg(coe2);
872 8365 : g2 = mkvec2(a, gsub(c, gaddgs(a, m))); /* 15.8.9 */
873 8365 : return FBaux2(v1, g1, coe1, m, z1, coe2, g2, v2, z1, Z, 1, prec);
874 : }
875 :
876 : /* 1 - z */
877 : static GEN
878 189 : F21taylor4(GEN a, GEN b, GEN c, GEN z, long prec)
879 : {
880 189 : GEN bma = gsub(c, gadd(a, b)), coe2, z1, z2, g1, g2, v1, v2, F = NULL;
881 : long m;
882 189 : z1 = gsubsg(1, z);
883 189 : if (!islong(bma,&m, prec))
884 : {
885 : GEN c1, a2, b2, c2;
886 119 : c1 = gsubsg(1, bma);
887 119 : a2 = gsub(c, a);
888 119 : b2 = gsub(c, b);
889 119 : c2 = gaddsg(1, bma);
890 119 : coe2 = gneg(gpow(z1, bma, prec));
891 119 : return FBaux1(mkvec3(a,b,c1), mkvec2(a2,b2), gen_1, mkvec3(a2,b2,c2),
892 : mkvec2(a,b), coe2, z1, bma, prec, prec);
893 : }
894 : /* c - (a+b) ~ m; if m < 0: 15.8.12 */
895 70 : if (m < 0) { m = -m; a = gsubgs(a, m); b = gsubgs(b, m); F = gpowgs(z1, -m); }
896 70 : v2 = g1 = mkvec2(gaddgs(a,m), gaddgs(b,m));
897 70 : v1 = g2 = mkvec2(a, b);
898 70 : z2 = gneg(z1); coe2 = gpowgs(z2, m); /* 15.8.10 */
899 70 : z1 = FBaux2(v1, g1, gen_1, m, z1, coe2, g2, v2, z1, z1, 4, prec);
900 70 : if (F) z1 = gmul(z1, F);
901 70 : return z1;
902 : }
903 :
904 : /* 1 - 1/z */
905 : static GEN
906 91 : F21taylor5(GEN a, GEN b, GEN c, GEN z, long prec)
907 : {
908 91 : GEN bma = gsub(c, gadd(a,b)), coe1, coe2, z1, g1, g2, v1, v2, z2, F = NULL;
909 : long m;
910 91 : if (!islong(bma,&m, prec))
911 : {
912 : GEN b1, c1, d1, e1, b2, c2;
913 49 : d1 = gsub(c, a);
914 49 : b1 = gsubsg(1, d1);
915 49 : c1 = gsubsg(1, bma);
916 49 : e1 = gsub(c, b);
917 49 : b2 = gsubsg(1, a);
918 49 : c2 = gaddsg(1, bma);
919 49 : coe1 = gpow(z, gneg(a), prec);
920 49 : coe2 = gneg(gmul(gpow(gsubsg(1, z), bma, prec), gpow(z, gneg(d1), prec)));
921 49 : z1 = gsubsg(1, ginv(z));
922 49 : return FBaux1(mkvec3(a,b1,c1), mkvec2(d1,e1), coe1,
923 : mkvec3(d1,b2,c2), mkvec2(a, b), coe2, z1, bma, prec, prec);
924 : }
925 : /* c - (a+b) ~ m; if m < 0: 15.8.12 */
926 42 : if (m < 0)
927 28 : { m = -m; a = gsubgs(a,m); b = gsubgs(b,m); F = gpowgs(gsubsg(1, z), -m); }
928 42 : g1 = mkvec2(gaddgs(a,m), gaddgs(b,m));
929 42 : v1 = mkvec2(a, gsubsg(1-m, b));
930 42 : v2 = mkvec2(gaddgs(a,m), gsubsg(1,b));
931 42 : z1 = gsubgs(ginv(z), 1);
932 42 : z2 = gneg(z1);
933 42 : g2 = mkvec2(a, b);
934 42 : coe1 = gpow(z, gneg(a), prec);
935 42 : coe2 = gmul(coe1, gpowgs(z2, m)); /* 15.8.11 */
936 42 : z1 = FBaux2(v1, g1, coe1, m, z1, coe2, g2, v2, z2, z1, 5, prec);
937 42 : if (F) z1 = gmul(z1, F);
938 42 : return z1;
939 : }
940 :
941 : /* 1 / z */
942 : static GEN
943 84 : F21taylor6(GEN a, GEN b, GEN c, GEN z, long prec)
944 : {
945 84 : GEN bma = gsub(b, a), cma, am, coe1, coe2, z1, g1, g2, v1, v2, z2, Z;
946 : long m;
947 84 : if (!islong(bma,&m, prec))
948 : {
949 : GEN e1, e2, b1, b2, c1, c2;
950 56 : b1 = gaddgs(gsub(a,c), 1);
951 56 : c1 = gsubsg(1, bma);
952 56 : e1 = gsub(c,a);
953 56 : b2 = gaddgs(gsub(b,c), 1);
954 56 : c2 = gaddgs(bma, 1);
955 56 : e2 = gsub(c, b);
956 56 : coe1 = gpow(gneg(z), gneg(a), prec);
957 56 : coe2 = gneg(gpow(gneg(z), gneg(b), prec));
958 56 : z1 = ginv(z);
959 56 : return FBaux1(mkvec3(a,b1,c1), mkvec2(b,e1), coe1,
960 : mkvec3(b,b2,c2), mkvec2(a,e2), coe2, z1, bma, prec, prec);
961 : }
962 : /* b - a ~ m */
963 28 : if (m < 0) { swap(a,b); m = -m; }
964 28 : cma = gsub(c, a); am = gaddgs(a,m);
965 28 : Z = gneg(z);
966 28 : coe1 = gpow(Z, gneg(a), prec);
967 28 : coe2 = gdiv(coe1, gpowgs(z, m));
968 28 : g1 = mkvec2(am, cma);
969 28 : v1 = mkvec2(a, gsubsg(1, cma));
970 28 : g2 = mkvec2(a, gsubgs(cma, m));
971 28 : v2 = mkvec2(am, gsubsg(m+1, cma));
972 28 : z2 = ginv(z);
973 28 : z1 = gneg(z2); /* 15.8.8 */
974 28 : return FBaux2(v1, g1, coe1, m, z1, coe2, g2, v2, z2, Z, 6, prec);
975 : }
976 :
977 : /* (new b, new c, new z): given by bind, cind, zind
978 : * case 1: (c-b,1+a-b,1/(1-z)) 15.8.3 15.8.6
979 : * case 2: (c-b,c,z/(z-1)) 15.8.1
980 : * case 3: (b,c,z)
981 : * case 4: (b,b-c+a+1,1-z) 15.8.4 15.8.7 15.8.10
982 : * case 5: (1+a-c,b-c+a+1,1-1/z) 15.8.5 15.8.7 15.8.11
983 : * case 6: (1+a-c,1+a-b,1/z) 15.8.2 15.8.6 */
984 : static GEN
985 24570 : F21taylorind(GEN a, GEN b, GEN c, GEN z, long ind, long prec)
986 : {
987 : GEN res;
988 24570 : switch (labs(ind))
989 : {
990 12572 : case 1: res = F21taylor1(a, b, c, z, prec); break;
991 11242 : case 2:
992 11242 : return gmul(Ftaylor(mkvec2(a, gsub(c, b)), mkvec(c),
993 : gdiv(z, gsubgs(z,1)), prec),
994 : gpow(gsubsg(1,z), gneg(a), prec));
995 392 : case 3:
996 392 : return Ftaylor(mkvec2(a,b), mkvec(c), z, prec);
997 189 : case 4: res = F21taylor4(a, b, c, z, prec); break;
998 91 : case 5: res = F21taylor5(a, b, c, z, prec); break;
999 84 : default:res = F21taylor6(a, b, c, z, prec); break;
1000 : }
1001 12936 : return gmul(ggamma(c, prec), res);
1002 : }
1003 :
1004 : static long
1005 3304 : hypersimplify(GEN *pn, GEN *pd)
1006 : {
1007 3304 : GEN n = *pn, d = *pd;
1008 3304 : long j, ld = lg(d), ln = lg(n);
1009 6300 : for (j = 1; j < ld; j++)
1010 : {
1011 3262 : GEN t = gel(d, j);
1012 : long k;
1013 8288 : for (k = 1; k < ln; k++)
1014 5292 : if (gequal(t, gel(n, k)))
1015 : {
1016 266 : *pd = vecsplice(d, j);
1017 266 : *pn = vecsplice(n, k); return hypersimplify(pd, pn) + 1;
1018 : }
1019 : }
1020 3038 : return 0;
1021 : }
1022 :
1023 : static GEN
1024 2086 : f_pochall(void *E, GEN n)
1025 : {
1026 2086 : GEN S, a, tmp, N = gel(E,1), D = gel(E,2);
1027 2086 : long j, prec = itou(gel(E,3));
1028 2086 : S = gen_0;
1029 9471 : for (j = 1; j < lg(N); j++)
1030 : {
1031 7385 : a = gel(N, j);
1032 7385 : tmp = gsub(glngamma(gadd(n, a), prec), glngamma(a, prec));
1033 7385 : S = gadd(S, tmp);
1034 : }
1035 7385 : for (j = 1; j < lg(D); j++)
1036 : {
1037 5299 : a = gel(D, j);
1038 5299 : tmp = gsub(glngamma(gadd(n, a), prec), glngamma(a, prec));
1039 5299 : S = gsub(S, tmp);
1040 : }
1041 2086 : return gexp(gsub(S, glngamma(gaddsg(1, n), prec)), prec);
1042 : }
1043 : static GEN
1044 371 : f_pochall_alt(void *E, GEN n)
1045 371 : { GEN z = f_pochall(E,n); return mpodd(n)? gneg(z): z; }
1046 : /* z = \pm1 */
1047 : static GEN
1048 63 : sumz(GEN N, GEN D, long z, long prec)
1049 : {
1050 63 : void *E = (void*)mkvec3(N, D, utoi(prec));
1051 : GEN tab, be;
1052 63 : if (z == -1) return sumalt(E, f_pochall_alt, gen_0, prec);
1053 56 : be = gsub(vecsum(D), vecsum(N)); check_hyp1(be);
1054 56 : tab = sumnummonieninit(be, NULL, gen_0, prec);
1055 56 : return sumnummonien(E, f_pochall, gen_0, tab, prec);
1056 : }
1057 :
1058 : static GEN
1059 6076 : hypergeom_arg(GEN x)
1060 : {
1061 6076 : if (!x) return cgetg(1,t_VEC);
1062 6048 : return (typ(x) == t_VEC)? x: mkvec(x);
1063 : }
1064 :
1065 : /* [x[i]+k, i=1..#v] */
1066 : static GEN
1067 6048 : RgV_z_add(GEN x, long k)
1068 : {
1069 6048 : if (!k) return x;
1070 13559 : pari_APPLY_same(gaddgs(gel(x,i), k));
1071 : }
1072 : static GEN
1073 161 : vp(GEN a, GEN p, GEN dft)
1074 : {
1075 161 : long v = gvaluation(a, p);
1076 161 : return v < 0? stoi(v): dft;
1077 : }
1078 : static GEN
1079 77 : Qp_hypergeom(GEN N, GEN D, GEN z)
1080 : {
1081 77 : pari_sp av = avma;
1082 77 : GEN r, S = gen_1, R = gen_1, p = padic_p(z), dft = ginv(subis(p, 1));
1083 77 : long l, i, prec = precp(z) + valp(z) + 1;
1084 :
1085 77 : r = gsub(stoi(valp(z)), dft);
1086 189 : l = lg(N); for (i = 1; i < l; i++) r = gadd(r, vp(gel(N,i), p, dft));
1087 126 : l = lg(D); for (i = 1; i < l; i++) r = gsub(r, vp(gel(D,i), p, dft));
1088 77 : if (gsigne(r) <= 0) pari_err(e_MISC, "divergent p-adic hypergeometric sum");
1089 56 : l = itou(gceil(gdivsg(prec, r)));
1090 2331 : for (i = 1; i <= l; i++)
1091 : {
1092 2275 : GEN u = vecprod(RgV_z_add(N, i-1)), v = vecprod(RgV_z_add(D, i-1));
1093 2275 : R = gmul(R, gmul(z, gdiv(u, gmulsg(i, v))));
1094 2275 : S = gadd(S, R);
1095 2275 : if (gc_needed(av,1))
1096 : {
1097 0 : if (DEBUGMEM>1) pari_warn(warnmem,"hypergeom, i = %ld / %ld", i,l);
1098 0 : (void)gc_all(av, 2, &R, &S);
1099 : }
1100 : }
1101 56 : return S;
1102 : }
1103 :
1104 : /* assume is_scalar_t(typ(z)) */
1105 : static GEN
1106 3325 : hypergeom_i(GEN N, GEN D, GEN z, long prec)
1107 : {
1108 : long nN, nD;
1109 3325 : if (typ(z) == t_PADIC) return Qp_hypergeom(N, D, z);
1110 3248 : if (gequal0(z)) return gen_1;
1111 3248 : nN = lg(N) - 1;
1112 3248 : nD = lg(D) - 1;
1113 3248 : if (nD >= (nN? nN: 2)) return Ftaylor(N, D, z, prec);
1114 2198 : if (nD == nN - 1 && nN >= 3)
1115 : {
1116 126 : GEN d = gsubsg(1, gabs(z,LOWDEFAULTPREC));
1117 126 : long ed = gexpo(d);
1118 : /* z in unit disc but "away" from unit circle */
1119 126 : if (gsigne(d) > 0 && ed > -prec2nbits(prec)/4
1120 21 : && (nN != 3 || ed > -15)) /* For 3F2 we can use integral */
1121 14 : return Ftaylor(N, D, z, prec);
1122 112 : if (gequal1(z)) return sumz(N, D, 1, prec);
1123 56 : if (gequalm1(z)) return sumz(N, D,-1, prec);
1124 : }
1125 2121 : switch (nN)
1126 : {
1127 112 : case 0:
1128 112 : if (nD == 0) return gexp(z, prec);
1129 21 : if (nD == 1) return F01(gel(D,1), z, prec);
1130 350 : case 1: return F10(gel(N, 1), z, prec);
1131 1547 : case 2:
1132 1547 : if (nD == 0) return F20(gel(N,1), gel(N,2), z, prec);
1133 1491 : if (nD == 1) return F21(gel(N,1), gel(N,2), gel(D,1), z, prec);
1134 : case 3:
1135 105 : if (nD == 0) break;
1136 105 : if (nD == 1) return F31(gel(N,1), gel(N,2), gel(N,3), gel(D,1), z, prec);
1137 49 : if (nD == 2) return F32(N, D, z, prec);
1138 : }
1139 7 : pari_err_IMPL("this hypergeometric function");
1140 : return NULL; /*LCOV_EXCL_LINE*/
1141 : }
1142 :
1143 : static GEN
1144 63 : serhypergeom(GEN N, GEN D, GEN y, long prec)
1145 : {
1146 63 : GEN Di, Ni, S = gen_1, R = gen_1, y0 = NULL;
1147 : long v, l, i;
1148 : pari_sp av;
1149 63 : if (!signe(y)) return gadd(gen_1, y);
1150 56 : v = valser(y); l = lg(y);
1151 56 : if (v < 0) pari_err_DOMAIN("hypergeom","valuation", "<", gen_0, y);
1152 56 : if (!v)
1153 : {
1154 28 : y0 = gel(y, 2);
1155 28 : if (!is_scalar_t(typ(y0))) pari_err_TYPE("hypergeom",y);
1156 28 : y = serchop0(y); l = 3 + (l - 3) / valser(y);
1157 28 : S = hypergeom(N, D, y0, prec);
1158 : }
1159 56 : av = avma; Ni = N; Di = D;
1160 805 : for (i = 1; i < l; i++)
1161 : {
1162 749 : R = gmul(R, gmul(y, gdiv(vecprod(Ni), gmulsg(i, vecprod(Di)))));
1163 : /* have to offset by 1 when y0 != NULL; keep Ni,Di to avoid recomputing
1164 : * in next loop */
1165 749 : Ni = RgV_z_add(N, i);
1166 749 : Di = RgV_z_add(D, i);
1167 749 : S = gadd(S, y0? gmul(R, hypergeom_i(Ni, Di, y0, prec)): R);
1168 749 : if (gc_needed(av,1))
1169 : {
1170 0 : if (DEBUGMEM>1) pari_warn(warnmem,"hypergeom, i = %ld / %ld", i,l-1);
1171 0 : (void)gc_all(av, 4, &S, &R, &Ni, &Di);
1172 : }
1173 : }
1174 56 : return S;
1175 : }
1176 :
1177 : GEN
1178 3038 : hypergeom(GEN N, GEN D, GEN y, long prec)
1179 : {
1180 3038 : pari_sp av = avma;
1181 : GEN z;
1182 : long j, n;
1183 3038 : N = hypergeom_arg(N);
1184 3038 : D = hypergeom_arg(D);
1185 3038 : n = hypersimplify(&N, &D);
1186 5866 : for (j = 1; j < lg(D); j++)
1187 2842 : if (isnegint(gel(D,j)))
1188 14 : pari_err_DOMAIN("hypergeom", stack_sprintf("b[%ld]", j + n),
1189 14 : "<=", gen_0, gel(D,j));
1190 3024 : if (is_scalar_t(typ(y)))
1191 2961 : return gc_GEN(av, hypergeom_i(N, D, y, prec));
1192 63 : if (!(z = toser_i(y))) pari_err_TYPE("hypergeom", y);
1193 63 : return gc_upto(av, serhypergeom(N, D, z, prec));
1194 : }
|
