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