Coverage for /usr/lib/python3/dist-packages/mpmath/functions/hypergeometric.py: 7%
905 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1from ..libmp.backend import xrange
2from .functions import defun, defun_wrapped
4def _check_need_perturb(ctx, terms, prec, discard_known_zeros):
5 perturb = recompute = False
6 extraprec = 0
7 discard = []
8 for term_index, term in enumerate(terms):
9 w_s, c_s, alpha_s, beta_s, a_s, b_s, z = term
10 have_singular_nongamma_weight = False
11 # Avoid division by zero in leading factors (TODO:
12 # also check for near division by zero?)
13 for k, w in enumerate(w_s):
14 if not w:
15 if ctx.re(c_s[k]) <= 0 and c_s[k]:
16 perturb = recompute = True
17 have_singular_nongamma_weight = True
18 pole_count = [0, 0, 0]
19 # Check for gamma and series poles and near-poles
20 for data_index, data in enumerate([alpha_s, beta_s, b_s]):
21 for i, x in enumerate(data):
22 n, d = ctx.nint_distance(x)
23 # Poles
24 if n > 0:
25 continue
26 if d == ctx.ninf:
27 # OK if we have a polynomial
28 # ------------------------------
29 ok = False
30 if data_index == 2:
31 for u in a_s:
32 if ctx.isnpint(u) and u >= int(n):
33 ok = True
34 break
35 if ok:
36 continue
37 pole_count[data_index] += 1
38 # ------------------------------
39 #perturb = recompute = True
40 #return perturb, recompute, extraprec
41 elif d < -4:
42 extraprec += -d
43 recompute = True
44 if discard_known_zeros and pole_count[1] > pole_count[0] + pole_count[2] \
45 and not have_singular_nongamma_weight:
46 discard.append(term_index)
47 elif sum(pole_count):
48 perturb = recompute = True
49 return perturb, recompute, extraprec, discard
51_hypercomb_msg = """
52hypercomb() failed to converge to the requested %i bits of accuracy
53using a working precision of %i bits. The function value may be zero or
54infinite; try passing zeroprec=N or infprec=M to bound finite values between
552^(-N) and 2^M. Otherwise try a higher maxprec or maxterms.
56"""
58@defun
59def hypercomb(ctx, function, params=[], discard_known_zeros=True, **kwargs):
60 orig = ctx.prec
61 sumvalue = ctx.zero
62 dist = ctx.nint_distance
63 ninf = ctx.ninf
64 orig_params = params[:]
65 verbose = kwargs.get('verbose', False)
66 maxprec = kwargs.get('maxprec', ctx._default_hyper_maxprec(orig))
67 kwargs['maxprec'] = maxprec # For calls to hypsum
68 zeroprec = kwargs.get('zeroprec')
69 infprec = kwargs.get('infprec')
70 perturbed_reference_value = None
71 hextra = 0
72 try:
73 while 1:
74 ctx.prec += 10
75 if ctx.prec > maxprec:
76 raise ValueError(_hypercomb_msg % (orig, ctx.prec))
77 orig2 = ctx.prec
78 params = orig_params[:]
79 terms = function(*params)
80 if verbose:
81 print()
82 print("ENTERING hypercomb main loop")
83 print("prec =", ctx.prec)
84 print("hextra", hextra)
85 perturb, recompute, extraprec, discard = \
86 _check_need_perturb(ctx, terms, orig, discard_known_zeros)
87 ctx.prec += extraprec
88 if perturb:
89 if "hmag" in kwargs:
90 hmag = kwargs["hmag"]
91 elif ctx._fixed_precision:
92 hmag = int(ctx.prec*0.3)
93 else:
94 hmag = orig + 10 + hextra
95 h = ctx.ldexp(ctx.one, -hmag)
96 ctx.prec = orig2 + 10 + hmag + 10
97 for k in range(len(params)):
98 params[k] += h
99 # Heuristically ensure that the perturbations
100 # are "independent" so that two perturbations
101 # don't accidentally cancel each other out
102 # in a subtraction.
103 h += h/(k+1)
104 if recompute:
105 terms = function(*params)
106 if discard_known_zeros:
107 terms = [term for (i, term) in enumerate(terms) if i not in discard]
108 if not terms:
109 return ctx.zero
110 evaluated_terms = []
111 for term_index, term_data in enumerate(terms):
112 w_s, c_s, alpha_s, beta_s, a_s, b_s, z = term_data
113 if verbose:
114 print()
115 print(" Evaluating term %i/%i : %iF%i" % \
116 (term_index+1, len(terms), len(a_s), len(b_s)))
117 print(" powers", ctx.nstr(w_s), ctx.nstr(c_s))
118 print(" gamma", ctx.nstr(alpha_s), ctx.nstr(beta_s))
119 print(" hyper", ctx.nstr(a_s), ctx.nstr(b_s))
120 print(" z", ctx.nstr(z))
121 #v = ctx.hyper(a_s, b_s, z, **kwargs)
122 #for a in alpha_s: v *= ctx.gamma(a)
123 #for b in beta_s: v *= ctx.rgamma(b)
124 #for w, c in zip(w_s, c_s): v *= ctx.power(w, c)
125 v = ctx.fprod([ctx.hyper(a_s, b_s, z, **kwargs)] + \
126 [ctx.gamma(a) for a in alpha_s] + \
127 [ctx.rgamma(b) for b in beta_s] + \
128 [ctx.power(w,c) for (w,c) in zip(w_s,c_s)])
129 if verbose:
130 print(" Value:", v)
131 evaluated_terms.append(v)
133 if len(terms) == 1 and (not perturb):
134 sumvalue = evaluated_terms[0]
135 break
137 if ctx._fixed_precision:
138 sumvalue = ctx.fsum(evaluated_terms)
139 break
141 sumvalue = ctx.fsum(evaluated_terms)
142 term_magnitudes = [ctx.mag(x) for x in evaluated_terms]
143 max_magnitude = max(term_magnitudes)
144 sum_magnitude = ctx.mag(sumvalue)
145 cancellation = max_magnitude - sum_magnitude
146 if verbose:
147 print()
148 print(" Cancellation:", cancellation, "bits")
149 print(" Increased precision:", ctx.prec - orig, "bits")
151 precision_ok = cancellation < ctx.prec - orig
153 if zeroprec is None:
154 zero_ok = False
155 else:
156 zero_ok = max_magnitude - ctx.prec < -zeroprec
157 if infprec is None:
158 inf_ok = False
159 else:
160 inf_ok = max_magnitude > infprec
162 if precision_ok and (not perturb) or ctx.isnan(cancellation):
163 break
164 elif precision_ok:
165 if perturbed_reference_value is None:
166 hextra += 20
167 perturbed_reference_value = sumvalue
168 continue
169 elif ctx.mag(sumvalue - perturbed_reference_value) <= \
170 ctx.mag(sumvalue) - orig:
171 break
172 elif zero_ok:
173 sumvalue = ctx.zero
174 break
175 elif inf_ok:
176 sumvalue = ctx.inf
177 break
178 elif 'hmag' in kwargs:
179 break
180 else:
181 hextra *= 2
182 perturbed_reference_value = sumvalue
183 # Increase precision
184 else:
185 increment = min(max(cancellation, orig//2), max(extraprec,orig))
186 ctx.prec += increment
187 if verbose:
188 print(" Must start over with increased precision")
189 continue
190 finally:
191 ctx.prec = orig
192 return +sumvalue
194@defun
195def hyper(ctx, a_s, b_s, z, **kwargs):
196 """
197 Hypergeometric function, general case.
198 """
199 z = ctx.convert(z)
200 p = len(a_s)
201 q = len(b_s)
202 a_s = [ctx._convert_param(a) for a in a_s]
203 b_s = [ctx._convert_param(b) for b in b_s]
204 # Reduce degree by eliminating common parameters
205 if kwargs.get('eliminate', True):
206 elim_nonpositive = kwargs.get('eliminate_all', False)
207 i = 0
208 while i < q and a_s:
209 b = b_s[i]
210 if b in a_s and (elim_nonpositive or not ctx.isnpint(b[0])):
211 a_s.remove(b)
212 b_s.remove(b)
213 p -= 1
214 q -= 1
215 else:
216 i += 1
217 # Handle special cases
218 if p == 0:
219 if q == 1: return ctx._hyp0f1(b_s, z, **kwargs)
220 elif q == 0: return ctx.exp(z)
221 elif p == 1:
222 if q == 1: return ctx._hyp1f1(a_s, b_s, z, **kwargs)
223 elif q == 2: return ctx._hyp1f2(a_s, b_s, z, **kwargs)
224 elif q == 0: return ctx._hyp1f0(a_s[0][0], z)
225 elif p == 2:
226 if q == 1: return ctx._hyp2f1(a_s, b_s, z, **kwargs)
227 elif q == 2: return ctx._hyp2f2(a_s, b_s, z, **kwargs)
228 elif q == 3: return ctx._hyp2f3(a_s, b_s, z, **kwargs)
229 elif q == 0: return ctx._hyp2f0(a_s, b_s, z, **kwargs)
230 elif p == q+1:
231 return ctx._hypq1fq(p, q, a_s, b_s, z, **kwargs)
232 elif p > q+1 and not kwargs.get('force_series'):
233 return ctx._hyp_borel(p, q, a_s, b_s, z, **kwargs)
234 coeffs, types = zip(*(a_s+b_s))
235 return ctx.hypsum(p, q, types, coeffs, z, **kwargs)
237@defun
238def hyp0f1(ctx,b,z,**kwargs):
239 return ctx.hyper([],[b],z,**kwargs)
241@defun
242def hyp1f1(ctx,a,b,z,**kwargs):
243 return ctx.hyper([a],[b],z,**kwargs)
245@defun
246def hyp1f2(ctx,a1,b1,b2,z,**kwargs):
247 return ctx.hyper([a1],[b1,b2],z,**kwargs)
249@defun
250def hyp2f1(ctx,a,b,c,z,**kwargs):
251 return ctx.hyper([a,b],[c],z,**kwargs)
253@defun
254def hyp2f2(ctx,a1,a2,b1,b2,z,**kwargs):
255 return ctx.hyper([a1,a2],[b1,b2],z,**kwargs)
257@defun
258def hyp2f3(ctx,a1,a2,b1,b2,b3,z,**kwargs):
259 return ctx.hyper([a1,a2],[b1,b2,b3],z,**kwargs)
261@defun
262def hyp2f0(ctx,a,b,z,**kwargs):
263 return ctx.hyper([a,b],[],z,**kwargs)
265@defun
266def hyp3f2(ctx,a1,a2,a3,b1,b2,z,**kwargs):
267 return ctx.hyper([a1,a2,a3],[b1,b2],z,**kwargs)
269@defun_wrapped
270def _hyp1f0(ctx, a, z):
271 return (1-z) ** (-a)
273@defun
274def _hyp0f1(ctx, b_s, z, **kwargs):
275 (b, btype), = b_s
276 if z:
277 magz = ctx.mag(z)
278 else:
279 magz = 0
280 if magz >= 8 and not kwargs.get('force_series'):
281 try:
282 # http://functions.wolfram.com/HypergeometricFunctions/
283 # Hypergeometric0F1/06/02/03/0004/
284 # TODO: handle the all-real case more efficiently!
285 # TODO: figure out how much precision is needed (exponential growth)
286 orig = ctx.prec
287 try:
288 ctx.prec += 12 + magz//2
289 def h():
290 w = ctx.sqrt(-z)
291 jw = ctx.j*w
292 u = 1/(4*jw)
293 c = ctx.mpq_1_2 - b
294 E = ctx.exp(2*jw)
295 T1 = ([-jw,E], [c,-1], [], [], [b-ctx.mpq_1_2, ctx.mpq_3_2-b], [], -u)
296 T2 = ([jw,E], [c,1], [], [], [b-ctx.mpq_1_2, ctx.mpq_3_2-b], [], u)
297 return T1, T2
298 v = ctx.hypercomb(h, [], force_series=True)
299 v = ctx.gamma(b)/(2*ctx.sqrt(ctx.pi))*v
300 finally:
301 ctx.prec = orig
302 if ctx._is_real_type(b) and ctx._is_real_type(z):
303 v = ctx._re(v)
304 return +v
305 except ctx.NoConvergence:
306 pass
307 return ctx.hypsum(0, 1, (btype,), [b], z, **kwargs)
309@defun
310def _hyp1f1(ctx, a_s, b_s, z, **kwargs):
311 (a, atype), = a_s
312 (b, btype), = b_s
313 if not z:
314 return ctx.one+z
315 magz = ctx.mag(z)
316 if magz >= 7 and not (ctx.isint(a) and ctx.re(a) <= 0):
317 if ctx.isinf(z):
318 if ctx.sign(a) == ctx.sign(b) == ctx.sign(z) == 1:
319 return ctx.inf
320 return ctx.nan * z
321 try:
322 try:
323 ctx.prec += magz
324 sector = ctx._im(z) < 0
325 def h(a,b):
326 if sector:
327 E = ctx.expjpi(ctx.fneg(a, exact=True))
328 else:
329 E = ctx.expjpi(a)
330 rz = 1/z
331 T1 = ([E,z], [1,-a], [b], [b-a], [a, 1+a-b], [], -rz)
332 T2 = ([ctx.exp(z),z], [1,a-b], [b], [a], [b-a, 1-a], [], rz)
333 return T1, T2
334 v = ctx.hypercomb(h, [a,b], force_series=True)
335 if ctx._is_real_type(a) and ctx._is_real_type(b) and ctx._is_real_type(z):
336 v = ctx._re(v)
337 return +v
338 except ctx.NoConvergence:
339 pass
340 finally:
341 ctx.prec -= magz
342 v = ctx.hypsum(1, 1, (atype, btype), [a, b], z, **kwargs)
343 return v
345def _hyp2f1_gosper(ctx,a,b,c,z,**kwargs):
346 # Use Gosper's recurrence
347 # See http://www.math.utexas.edu/pipermail/maxima/2006/000126.html
348 _a,_b,_c,_z = a, b, c, z
349 orig = ctx.prec
350 maxprec = kwargs.get('maxprec', 100*orig)
351 extra = 10
352 while 1:
353 ctx.prec = orig + extra
354 #a = ctx.convert(_a)
355 #b = ctx.convert(_b)
356 #c = ctx.convert(_c)
357 z = ctx.convert(_z)
358 d = ctx.mpf(0)
359 e = ctx.mpf(1)
360 f = ctx.mpf(0)
361 k = 0
362 # Common subexpression elimination, unfortunately making
363 # things a bit unreadable. The formula is quite messy to begin
364 # with, though...
365 abz = a*b*z
366 ch = c * ctx.mpq_1_2
367 c1h = (c+1) * ctx.mpq_1_2
368 nz = 1-z
369 g = z/nz
370 abg = a*b*g
371 cba = c-b-a
372 z2 = z-2
373 tol = -ctx.prec - 10
374 nstr = ctx.nstr
375 nprint = ctx.nprint
376 mag = ctx.mag
377 maxmag = ctx.ninf
378 while 1:
379 kch = k+ch
380 kakbz = (k+a)*(k+b)*z / (4*(k+1)*kch*(k+c1h))
381 d1 = kakbz*(e-(k+cba)*d*g)
382 e1 = kakbz*(d*abg+(k+c)*e)
383 ft = d*(k*(cba*z+k*z2-c)-abz)/(2*kch*nz)
384 f1 = f + e - ft
385 maxmag = max(maxmag, mag(f1))
386 if mag(f1-f) < tol:
387 break
388 d, e, f = d1, e1, f1
389 k += 1
390 cancellation = maxmag - mag(f1)
391 if cancellation < extra:
392 break
393 else:
394 extra += cancellation
395 if extra > maxprec:
396 raise ctx.NoConvergence
397 return f1
399@defun
400def _hyp2f1(ctx, a_s, b_s, z, **kwargs):
401 (a, atype), (b, btype) = a_s
402 (c, ctype), = b_s
403 if z == 1:
404 # TODO: the following logic can be simplified
405 convergent = ctx.re(c-a-b) > 0
406 finite = (ctx.isint(a) and a <= 0) or (ctx.isint(b) and b <= 0)
407 zerodiv = ctx.isint(c) and c <= 0 and not \
408 ((ctx.isint(a) and c <= a <= 0) or (ctx.isint(b) and c <= b <= 0))
409 #print "bz", a, b, c, z, convergent, finite, zerodiv
410 # Gauss's theorem gives the value if convergent
411 if (convergent or finite) and not zerodiv:
412 return ctx.gammaprod([c, c-a-b], [c-a, c-b], _infsign=True)
413 # Otherwise, there is a pole and we take the
414 # sign to be that when approaching from below
415 # XXX: this evaluation is not necessarily correct in all cases
416 return ctx.hyp2f1(a,b,c,1-ctx.eps*2) * ctx.inf
418 # Equal to 1 (first term), unless there is a subsequent
419 # division by zero
420 if not z:
421 # Division by zero but power of z is higher than
422 # first order so cancels
423 if c or a == 0 or b == 0:
424 return 1+z
425 # Indeterminate
426 return ctx.nan
428 # Hit zero denominator unless numerator goes to 0 first
429 if ctx.isint(c) and c <= 0:
430 if (ctx.isint(a) and c <= a <= 0) or \
431 (ctx.isint(b) and c <= b <= 0):
432 pass
433 else:
434 # Pole in series
435 return ctx.inf
437 absz = abs(z)
439 # Fast case: standard series converges rapidly,
440 # possibly in finitely many terms
441 if absz <= 0.8 or (ctx.isint(a) and a <= 0 and a >= -1000) or \
442 (ctx.isint(b) and b <= 0 and b >= -1000):
443 return ctx.hypsum(2, 1, (atype, btype, ctype), [a, b, c], z, **kwargs)
445 orig = ctx.prec
446 try:
447 ctx.prec += 10
449 # Use 1/z transformation
450 if absz >= 1.3:
451 def h(a,b):
452 t = ctx.mpq_1-c; ab = a-b; rz = 1/z
453 T1 = ([-z],[-a], [c,-ab],[b,c-a], [a,t+a],[ctx.mpq_1+ab], rz)
454 T2 = ([-z],[-b], [c,ab],[a,c-b], [b,t+b],[ctx.mpq_1-ab], rz)
455 return T1, T2
456 v = ctx.hypercomb(h, [a,b], **kwargs)
458 # Use 1-z transformation
459 elif abs(1-z) <= 0.75:
460 def h(a,b):
461 t = c-a-b; ca = c-a; cb = c-b; rz = 1-z
462 T1 = [], [], [c,t], [ca,cb], [a,b], [1-t], rz
463 T2 = [rz], [t], [c,a+b-c], [a,b], [ca,cb], [1+t], rz
464 return T1, T2
465 v = ctx.hypercomb(h, [a,b], **kwargs)
467 # Use z/(z-1) transformation
468 elif abs(z/(z-1)) <= 0.75:
469 v = ctx.hyp2f1(a, c-b, c, z/(z-1)) / (1-z)**a
471 # Remaining part of unit circle
472 else:
473 v = _hyp2f1_gosper(ctx,a,b,c,z,**kwargs)
475 finally:
476 ctx.prec = orig
477 return +v
479@defun
480def _hypq1fq(ctx, p, q, a_s, b_s, z, **kwargs):
481 r"""
482 Evaluates 3F2, 4F3, 5F4, ...
483 """
484 a_s, a_types = zip(*a_s)
485 b_s, b_types = zip(*b_s)
486 a_s = list(a_s)
487 b_s = list(b_s)
488 absz = abs(z)
489 ispoly = False
490 for a in a_s:
491 if ctx.isint(a) and a <= 0:
492 ispoly = True
493 break
494 # Direct summation
495 if absz < 1 or ispoly:
496 try:
497 return ctx.hypsum(p, q, a_types+b_types, a_s+b_s, z, **kwargs)
498 except ctx.NoConvergence:
499 if absz > 1.1 or ispoly:
500 raise
501 # Use expansion at |z-1| -> 0.
502 # Reference: Wolfgang Buhring, "Generalized Hypergeometric Functions at
503 # Unit Argument", Proc. Amer. Math. Soc., Vol. 114, No. 1 (Jan. 1992),
504 # pp.145-153
505 # The current implementation has several problems:
506 # 1. We only implement it for 3F2. The expansion coefficients are
507 # given by extremely messy nested sums in the higher degree cases
508 # (see reference). Is efficient sequential generation of the coefficients
509 # possible in the > 3F2 case?
510 # 2. Although the series converges, it may do so slowly, so we need
511 # convergence acceleration. The acceleration implemented by
512 # nsum does not always help, so results returned are sometimes
513 # inaccurate! Can we do better?
514 # 3. We should check conditions for convergence, and possibly
515 # do a better job of cancelling out gamma poles if possible.
516 if z == 1:
517 # XXX: should also check for division by zero in the
518 # denominator of the series (cf. hyp2f1)
519 S = ctx.re(sum(b_s)-sum(a_s))
520 if S <= 0:
521 #return ctx.hyper(a_s, b_s, 1-ctx.eps*2, **kwargs) * ctx.inf
522 return ctx.hyper(a_s, b_s, 0.9, **kwargs) * ctx.inf
523 if (p,q) == (3,2) and abs(z-1) < 0.05: # and kwargs.get('sum1')
524 #print "Using alternate summation (experimental)"
525 a1,a2,a3 = a_s
526 b1,b2 = b_s
527 u = b1+b2-a3
528 initial = ctx.gammaprod([b2-a3,b1-a3,a1,a2],[b2-a3,b1-a3,1,u])
529 def term(k, _cache={0:initial}):
530 u = b1+b2-a3+k
531 if k in _cache:
532 t = _cache[k]
533 else:
534 t = _cache[k-1]
535 t *= (b1+k-a3-1)*(b2+k-a3-1)
536 t /= k*(u-1)
537 _cache[k] = t
538 return t * ctx.hyp2f1(a1,a2,u,z)
539 try:
540 S = ctx.nsum(term, [0,ctx.inf], verbose=kwargs.get('verbose'),
541 strict=kwargs.get('strict', True))
542 return S * ctx.gammaprod([b1,b2],[a1,a2,a3])
543 except ctx.NoConvergence:
544 pass
545 # Try to use convergence acceleration on and close to the unit circle.
546 # Problem: the convergence acceleration degenerates as |z-1| -> 0,
547 # except for special cases. Everywhere else, the Shanks transformation
548 # is very efficient.
549 if absz < 1.1 and ctx._re(z) <= 1:
551 def term(kk, _cache={0:ctx.one}):
552 k = int(kk)
553 if k != kk:
554 t = z ** ctx.mpf(kk) / ctx.fac(kk)
555 for a in a_s: t *= ctx.rf(a,kk)
556 for b in b_s: t /= ctx.rf(b,kk)
557 return t
558 if k in _cache:
559 return _cache[k]
560 t = term(k-1)
561 m = k-1
562 for j in xrange(p): t *= (a_s[j]+m)
563 for j in xrange(q): t /= (b_s[j]+m)
564 t *= z
565 t /= k
566 _cache[k] = t
567 return t
569 sum_method = kwargs.get('sum_method', 'r+s+e')
571 try:
572 return ctx.nsum(term, [0,ctx.inf], verbose=kwargs.get('verbose'),
573 strict=kwargs.get('strict', True),
574 method=sum_method.replace('e',''))
575 except ctx.NoConvergence:
576 if 'e' not in sum_method:
577 raise
578 pass
580 if kwargs.get('verbose'):
581 print("Attempting Euler-Maclaurin summation")
584 """
585 Somewhat slower version (one diffs_exp for each factor).
586 However, this would be faster with fast direct derivatives
587 of the gamma function.
589 def power_diffs(k0):
590 r = 0
591 l = ctx.log(z)
592 while 1:
593 yield z**ctx.mpf(k0) * l**r
594 r += 1
596 def loggamma_diffs(x, reciprocal=False):
597 sign = (-1) ** reciprocal
598 yield sign * ctx.loggamma(x)
599 i = 0
600 while 1:
601 yield sign * ctx.psi(i,x)
602 i += 1
604 def hyper_diffs(k0):
605 b2 = b_s + [1]
606 A = [ctx.diffs_exp(loggamma_diffs(a+k0)) for a in a_s]
607 B = [ctx.diffs_exp(loggamma_diffs(b+k0,True)) for b in b2]
608 Z = [power_diffs(k0)]
609 C = ctx.gammaprod([b for b in b2], [a for a in a_s])
610 for d in ctx.diffs_prod(A + B + Z):
611 v = C * d
612 yield v
613 """
615 def log_diffs(k0):
616 b2 = b_s + [1]
617 yield sum(ctx.loggamma(a+k0) for a in a_s) - \
618 sum(ctx.loggamma(b+k0) for b in b2) + k0*ctx.log(z)
619 i = 0
620 while 1:
621 v = sum(ctx.psi(i,a+k0) for a in a_s) - \
622 sum(ctx.psi(i,b+k0) for b in b2)
623 if i == 0:
624 v += ctx.log(z)
625 yield v
626 i += 1
628 def hyper_diffs(k0):
629 C = ctx.gammaprod([b for b in b_s], [a for a in a_s])
630 for d in ctx.diffs_exp(log_diffs(k0)):
631 v = C * d
632 yield v
634 tol = ctx.eps / 1024
635 prec = ctx.prec
636 try:
637 trunc = 50 * ctx.dps
638 ctx.prec += 20
639 for i in xrange(5):
640 head = ctx.fsum(term(k) for k in xrange(trunc))
641 tail, err = ctx.sumem(term, [trunc, ctx.inf], tol=tol,
642 adiffs=hyper_diffs(trunc),
643 verbose=kwargs.get('verbose'),
644 error=True,
645 _fast_abort=True)
646 if err < tol:
647 v = head + tail
648 break
649 trunc *= 2
650 # Need to increase precision because calculation of
651 # derivatives may be inaccurate
652 ctx.prec += ctx.prec//2
653 if i == 4:
654 raise ctx.NoConvergence(\
655 "Euler-Maclaurin summation did not converge")
656 finally:
657 ctx.prec = prec
658 return +v
660 # Use 1/z transformation
661 # http://functions.wolfram.com/HypergeometricFunctions/
662 # HypergeometricPFQ/06/01/05/02/0004/
663 def h(*args):
664 a_s = list(args[:p])
665 b_s = list(args[p:])
666 Ts = []
667 recz = ctx.one/z
668 negz = ctx.fneg(z, exact=True)
669 for k in range(q+1):
670 ak = a_s[k]
671 C = [negz]
672 Cp = [-ak]
673 Gn = b_s + [ak] + [a_s[j]-ak for j in range(q+1) if j != k]
674 Gd = a_s + [b_s[j]-ak for j in range(q)]
675 Fn = [ak] + [ak-b_s[j]+1 for j in range(q)]
676 Fd = [1-a_s[j]+ak for j in range(q+1) if j != k]
677 Ts.append((C, Cp, Gn, Gd, Fn, Fd, recz))
678 return Ts
679 return ctx.hypercomb(h, a_s+b_s, **kwargs)
681@defun
682def _hyp_borel(ctx, p, q, a_s, b_s, z, **kwargs):
683 if a_s:
684 a_s, a_types = zip(*a_s)
685 a_s = list(a_s)
686 else:
687 a_s, a_types = [], ()
688 if b_s:
689 b_s, b_types = zip(*b_s)
690 b_s = list(b_s)
691 else:
692 b_s, b_types = [], ()
693 kwargs['maxterms'] = kwargs.get('maxterms', ctx.prec)
694 try:
695 return ctx.hypsum(p, q, a_types+b_types, a_s+b_s, z, **kwargs)
696 except ctx.NoConvergence:
697 pass
698 prec = ctx.prec
699 try:
700 tol = kwargs.get('asymp_tol', ctx.eps/4)
701 ctx.prec += 10
702 # hypsum is has a conservative tolerance. So we try again:
703 def term(k, cache={0:ctx.one}):
704 if k in cache:
705 return cache[k]
706 t = term(k-1)
707 for a in a_s: t *= (a+(k-1))
708 for b in b_s: t /= (b+(k-1))
709 t *= z
710 t /= k
711 cache[k] = t
712 return t
713 s = ctx.one
714 for k in xrange(1, ctx.prec):
715 t = term(k)
716 s += t
717 if abs(t) <= tol:
718 return s
719 finally:
720 ctx.prec = prec
721 if p <= q+3:
722 contour = kwargs.get('contour')
723 if not contour:
724 if ctx.arg(z) < 0.25:
725 u = z / max(1, abs(z))
726 if ctx.arg(z) >= 0:
727 contour = [0, 2j, (2j+2)/u, 2/u, ctx.inf]
728 else:
729 contour = [0, -2j, (-2j+2)/u, 2/u, ctx.inf]
730 #contour = [0, 2j/z, 2/z, ctx.inf]
731 #contour = [0, 2j, 2/z, ctx.inf]
732 #contour = [0, 2j, ctx.inf]
733 else:
734 contour = [0, ctx.inf]
735 quad_kwargs = kwargs.get('quad_kwargs', {})
736 def g(t):
737 return ctx.exp(-t)*ctx.hyper(a_s, b_s+[1], t*z)
738 I, err = ctx.quad(g, contour, error=True, **quad_kwargs)
739 if err <= abs(I)*ctx.eps*8:
740 return I
741 raise ctx.NoConvergence
744@defun
745def _hyp2f2(ctx, a_s, b_s, z, **kwargs):
746 (a1, a1type), (a2, a2type) = a_s
747 (b1, b1type), (b2, b2type) = b_s
749 absz = abs(z)
750 magz = ctx.mag(z)
751 orig = ctx.prec
753 # Asymptotic expansion is ~ exp(z)
754 asymp_extraprec = magz
756 # Asymptotic series is in terms of 3F1
757 can_use_asymptotic = (not kwargs.get('force_series')) and \
758 (ctx.mag(absz) > 3)
760 # TODO: much of the following could be shared with 2F3 instead of
761 # copypasted
762 if can_use_asymptotic:
763 #print "using asymp"
764 try:
765 try:
766 ctx.prec += asymp_extraprec
767 # http://functions.wolfram.com/HypergeometricFunctions/
768 # Hypergeometric2F2/06/02/02/0002/
769 def h(a1,a2,b1,b2):
770 X = a1+a2-b1-b2
771 A2 = a1+a2
772 B2 = b1+b2
773 c = {}
774 c[0] = ctx.one
775 c[1] = (A2-1)*X+b1*b2-a1*a2
776 s1 = 0
777 k = 0
778 tprev = 0
779 while 1:
780 if k not in c:
781 uu1 = 1-B2+2*a1+a1**2+2*a2+a2**2-A2*B2+a1*a2+b1*b2+(2*B2-3*(A2+1))*k+2*k**2
782 uu2 = (k-A2+b1-1)*(k-A2+b2-1)*(k-X-2)
783 c[k] = ctx.one/k * (uu1*c[k-1]-uu2*c[k-2])
784 t1 = c[k] * z**(-k)
785 if abs(t1) < 0.1*ctx.eps:
786 #print "Convergence :)"
787 break
788 # Quit if the series doesn't converge quickly enough
789 if k > 5 and abs(tprev) / abs(t1) < 1.5:
790 #print "No convergence :("
791 raise ctx.NoConvergence
792 s1 += t1
793 tprev = t1
794 k += 1
795 S = ctx.exp(z)*s1
796 T1 = [z,S], [X,1], [b1,b2],[a1,a2],[],[],0
797 T2 = [-z],[-a1],[b1,b2,a2-a1],[a2,b1-a1,b2-a1],[a1,a1-b1+1,a1-b2+1],[a1-a2+1],-1/z
798 T3 = [-z],[-a2],[b1,b2,a1-a2],[a1,b1-a2,b2-a2],[a2,a2-b1+1,a2-b2+1],[-a1+a2+1],-1/z
799 return T1, T2, T3
800 v = ctx.hypercomb(h, [a1,a2,b1,b2], force_series=True, maxterms=4*ctx.prec)
801 if sum(ctx._is_real_type(u) for u in [a1,a2,b1,b2,z]) == 5:
802 v = ctx.re(v)
803 return v
804 except ctx.NoConvergence:
805 pass
806 finally:
807 ctx.prec = orig
809 return ctx.hypsum(2, 2, (a1type, a2type, b1type, b2type), [a1, a2, b1, b2], z, **kwargs)
813@defun
814def _hyp1f2(ctx, a_s, b_s, z, **kwargs):
815 (a1, a1type), = a_s
816 (b1, b1type), (b2, b2type) = b_s
818 absz = abs(z)
819 magz = ctx.mag(z)
820 orig = ctx.prec
822 # Asymptotic expansion is ~ exp(sqrt(z))
823 asymp_extraprec = z and magz//2
825 # Asymptotic series is in terms of 3F0
826 can_use_asymptotic = (not kwargs.get('force_series')) and \
827 (ctx.mag(absz) > 19) and \
828 (ctx.sqrt(absz) > 1.5*orig) # and \
829 # ctx._hyp_check_convergence([a1, a1-b1+1, a1-b2+1], [],
830 # 1/absz, orig+40+asymp_extraprec)
832 # TODO: much of the following could be shared with 2F3 instead of
833 # copypasted
834 if can_use_asymptotic:
835 #print "using asymp"
836 try:
837 try:
838 ctx.prec += asymp_extraprec
839 # http://functions.wolfram.com/HypergeometricFunctions/
840 # Hypergeometric1F2/06/02/03/
841 def h(a1,b1,b2):
842 X = ctx.mpq_1_2*(a1-b1-b2+ctx.mpq_1_2)
843 c = {}
844 c[0] = ctx.one
845 c[1] = 2*(ctx.mpq_1_4*(3*a1+b1+b2-2)*(a1-b1-b2)+b1*b2-ctx.mpq_3_16)
846 c[2] = 2*(b1*b2+ctx.mpq_1_4*(a1-b1-b2)*(3*a1+b1+b2-2)-ctx.mpq_3_16)**2+\
847 ctx.mpq_1_16*(-16*(2*a1-3)*b1*b2 + \
848 4*(a1-b1-b2)*(-8*a1**2+11*a1+b1+b2-2)-3)
849 s1 = 0
850 s2 = 0
851 k = 0
852 tprev = 0
853 while 1:
854 if k not in c:
855 uu1 = (3*k**2+(-6*a1+2*b1+2*b2-4)*k + 3*a1**2 - \
856 (b1-b2)**2 - 2*a1*(b1+b2-2) + ctx.mpq_1_4)
857 uu2 = (k-a1+b1-b2-ctx.mpq_1_2)*(k-a1-b1+b2-ctx.mpq_1_2)*\
858 (k-a1+b1+b2-ctx.mpq_5_2)
859 c[k] = ctx.one/(2*k)*(uu1*c[k-1]-uu2*c[k-2])
860 w = c[k] * (-z)**(-0.5*k)
861 t1 = (-ctx.j)**k * ctx.mpf(2)**(-k) * w
862 t2 = ctx.j**k * ctx.mpf(2)**(-k) * w
863 if abs(t1) < 0.1*ctx.eps:
864 #print "Convergence :)"
865 break
866 # Quit if the series doesn't converge quickly enough
867 if k > 5 and abs(tprev) / abs(t1) < 1.5:
868 #print "No convergence :("
869 raise ctx.NoConvergence
870 s1 += t1
871 s2 += t2
872 tprev = t1
873 k += 1
874 S = ctx.expj(ctx.pi*X+2*ctx.sqrt(-z))*s1 + \
875 ctx.expj(-(ctx.pi*X+2*ctx.sqrt(-z)))*s2
876 T1 = [0.5*S, ctx.pi, -z], [1, -0.5, X], [b1, b2], [a1],\
877 [], [], 0
878 T2 = [-z], [-a1], [b1,b2],[b1-a1,b2-a1], \
879 [a1,a1-b1+1,a1-b2+1], [], 1/z
880 return T1, T2
881 v = ctx.hypercomb(h, [a1,b1,b2], force_series=True, maxterms=4*ctx.prec)
882 if sum(ctx._is_real_type(u) for u in [a1,b1,b2,z]) == 4:
883 v = ctx.re(v)
884 return v
885 except ctx.NoConvergence:
886 pass
887 finally:
888 ctx.prec = orig
890 #print "not using asymp"
891 return ctx.hypsum(1, 2, (a1type, b1type, b2type), [a1, b1, b2], z, **kwargs)
895@defun
896def _hyp2f3(ctx, a_s, b_s, z, **kwargs):
897 (a1, a1type), (a2, a2type) = a_s
898 (b1, b1type), (b2, b2type), (b3, b3type) = b_s
900 absz = abs(z)
901 magz = ctx.mag(z)
903 # Asymptotic expansion is ~ exp(sqrt(z))
904 asymp_extraprec = z and magz//2
905 orig = ctx.prec
907 # Asymptotic series is in terms of 4F1
908 # The square root below empirically provides a plausible criterion
909 # for the leading series to converge
910 can_use_asymptotic = (not kwargs.get('force_series')) and \
911 (ctx.mag(absz) > 19) and (ctx.sqrt(absz) > 1.5*orig)
913 if can_use_asymptotic:
914 #print "using asymp"
915 try:
916 try:
917 ctx.prec += asymp_extraprec
918 # http://functions.wolfram.com/HypergeometricFunctions/
919 # Hypergeometric2F3/06/02/03/01/0002/
920 def h(a1,a2,b1,b2,b3):
921 X = ctx.mpq_1_2*(a1+a2-b1-b2-b3+ctx.mpq_1_2)
922 A2 = a1+a2
923 B3 = b1+b2+b3
924 A = a1*a2
925 B = b1*b2+b3*b2+b1*b3
926 R = b1*b2*b3
927 c = {}
928 c[0] = ctx.one
929 c[1] = 2*(B - A + ctx.mpq_1_4*(3*A2+B3-2)*(A2-B3) - ctx.mpq_3_16)
930 c[2] = ctx.mpq_1_2*c[1]**2 + ctx.mpq_1_16*(-16*(2*A2-3)*(B-A) + 32*R +\
931 4*(-8*A2**2 + 11*A2 + 8*A + B3 - 2)*(A2-B3)-3)
932 s1 = 0
933 s2 = 0
934 k = 0
935 tprev = 0
936 while 1:
937 if k not in c:
938 uu1 = (k-2*X-3)*(k-2*X-2*b1-1)*(k-2*X-2*b2-1)*\
939 (k-2*X-2*b3-1)
940 uu2 = (4*(k-1)**3 - 6*(4*X+B3)*(k-1)**2 + \
941 2*(24*X**2+12*B3*X+4*B+B3-1)*(k-1) - 32*X**3 - \
942 24*B3*X**2 - 4*B - 8*R - 4*(4*B+B3-1)*X + 2*B3-1)
943 uu3 = (5*(k-1)**2+2*(-10*X+A2-3*B3+3)*(k-1)+2*c[1])
944 c[k] = ctx.one/(2*k)*(uu1*c[k-3]-uu2*c[k-2]+uu3*c[k-1])
945 w = c[k] * ctx.power(-z, -0.5*k)
946 t1 = (-ctx.j)**k * ctx.mpf(2)**(-k) * w
947 t2 = ctx.j**k * ctx.mpf(2)**(-k) * w
948 if abs(t1) < 0.1*ctx.eps:
949 break
950 # Quit if the series doesn't converge quickly enough
951 if k > 5 and abs(tprev) / abs(t1) < 1.5:
952 raise ctx.NoConvergence
953 s1 += t1
954 s2 += t2
955 tprev = t1
956 k += 1
957 S = ctx.expj(ctx.pi*X+2*ctx.sqrt(-z))*s1 + \
958 ctx.expj(-(ctx.pi*X+2*ctx.sqrt(-z)))*s2
959 T1 = [0.5*S, ctx.pi, -z], [1, -0.5, X], [b1, b2, b3], [a1, a2],\
960 [], [], 0
961 T2 = [-z], [-a1], [b1,b2,b3,a2-a1],[a2,b1-a1,b2-a1,b3-a1], \
962 [a1,a1-b1+1,a1-b2+1,a1-b3+1], [a1-a2+1], 1/z
963 T3 = [-z], [-a2], [b1,b2,b3,a1-a2],[a1,b1-a2,b2-a2,b3-a2], \
964 [a2,a2-b1+1,a2-b2+1,a2-b3+1],[-a1+a2+1], 1/z
965 return T1, T2, T3
966 v = ctx.hypercomb(h, [a1,a2,b1,b2,b3], force_series=True, maxterms=4*ctx.prec)
967 if sum(ctx._is_real_type(u) for u in [a1,a2,b1,b2,b3,z]) == 6:
968 v = ctx.re(v)
969 return v
970 except ctx.NoConvergence:
971 pass
972 finally:
973 ctx.prec = orig
975 return ctx.hypsum(2, 3, (a1type, a2type, b1type, b2type, b3type), [a1, a2, b1, b2, b3], z, **kwargs)
977@defun
978def _hyp2f0(ctx, a_s, b_s, z, **kwargs):
979 (a, atype), (b, btype) = a_s
980 # We want to try aggressively to use the asymptotic expansion,
981 # and fall back only when absolutely necessary
982 try:
983 kwargsb = kwargs.copy()
984 kwargsb['maxterms'] = kwargsb.get('maxterms', ctx.prec)
985 return ctx.hypsum(2, 0, (atype,btype), [a,b], z, **kwargsb)
986 except ctx.NoConvergence:
987 if kwargs.get('force_series'):
988 raise
989 pass
990 def h(a, b):
991 w = ctx.sinpi(b)
992 rz = -1/z
993 T1 = ([ctx.pi,w,rz],[1,-1,a],[],[a-b+1,b],[a],[b],rz)
994 T2 = ([-ctx.pi,w,rz],[1,-1,1+a-b],[],[a,2-b],[a-b+1],[2-b],rz)
995 return T1, T2
996 return ctx.hypercomb(h, [a, 1+a-b], **kwargs)
998@defun
999def meijerg(ctx, a_s, b_s, z, r=1, series=None, **kwargs):
1000 an, ap = a_s
1001 bm, bq = b_s
1002 n = len(an)
1003 p = n + len(ap)
1004 m = len(bm)
1005 q = m + len(bq)
1006 a = an+ap
1007 b = bm+bq
1008 a = [ctx.convert(_) for _ in a]
1009 b = [ctx.convert(_) for _ in b]
1010 z = ctx.convert(z)
1011 if series is None:
1012 if p < q: series = 1
1013 if p > q: series = 2
1014 if p == q:
1015 if m+n == p and abs(z) > 1:
1016 series = 2
1017 else:
1018 series = 1
1019 if kwargs.get('verbose'):
1020 print("Meijer G m,n,p,q,series =", m,n,p,q,series)
1021 if series == 1:
1022 def h(*args):
1023 a = args[:p]
1024 b = args[p:]
1025 terms = []
1026 for k in range(m):
1027 bases = [z]
1028 expts = [b[k]/r]
1029 gn = [b[j]-b[k] for j in range(m) if j != k]
1030 gn += [1-a[j]+b[k] for j in range(n)]
1031 gd = [a[j]-b[k] for j in range(n,p)]
1032 gd += [1-b[j]+b[k] for j in range(m,q)]
1033 hn = [1-a[j]+b[k] for j in range(p)]
1034 hd = [1-b[j]+b[k] for j in range(q) if j != k]
1035 hz = (-ctx.one)**(p-m-n) * z**(ctx.one/r)
1036 terms.append((bases, expts, gn, gd, hn, hd, hz))
1037 return terms
1038 else:
1039 def h(*args):
1040 a = args[:p]
1041 b = args[p:]
1042 terms = []
1043 for k in range(n):
1044 bases = [z]
1045 if r == 1:
1046 expts = [a[k]-1]
1047 else:
1048 expts = [(a[k]-1)/ctx.convert(r)]
1049 gn = [a[k]-a[j] for j in range(n) if j != k]
1050 gn += [1-a[k]+b[j] for j in range(m)]
1051 gd = [a[k]-b[j] for j in range(m,q)]
1052 gd += [1-a[k]+a[j] for j in range(n,p)]
1053 hn = [1-a[k]+b[j] for j in range(q)]
1054 hd = [1+a[j]-a[k] for j in range(p) if j != k]
1055 hz = (-ctx.one)**(q-m-n) / z**(ctx.one/r)
1056 terms.append((bases, expts, gn, gd, hn, hd, hz))
1057 return terms
1058 return ctx.hypercomb(h, a+b, **kwargs)
1060@defun_wrapped
1061def appellf1(ctx,a,b1,b2,c,x,y,**kwargs):
1062 # Assume x smaller
1063 # We will use x for the outer loop
1064 if abs(x) > abs(y):
1065 x, y = y, x
1066 b1, b2 = b2, b1
1067 def ok(x):
1068 return abs(x) < 0.99
1069 # Finite cases
1070 if ctx.isnpint(a):
1071 pass
1072 elif ctx.isnpint(b1):
1073 pass
1074 elif ctx.isnpint(b2):
1075 x, y, b1, b2 = y, x, b2, b1
1076 else:
1077 #print x, y
1078 # Note: ok if |y| > 1, because
1079 # 2F1 implements analytic continuation
1080 if not ok(x):
1081 u1 = (x-y)/(x-1)
1082 if not ok(u1):
1083 raise ValueError("Analytic continuation not implemented")
1084 #print "Using analytic continuation"
1085 return (1-x)**(-b1)*(1-y)**(c-a-b2)*\
1086 ctx.appellf1(c-a,b1,c-b1-b2,c,u1,y,**kwargs)
1087 return ctx.hyper2d({'m+n':[a],'m':[b1],'n':[b2]}, {'m+n':[c]}, x,y, **kwargs)
1089@defun
1090def appellf2(ctx,a,b1,b2,c1,c2,x,y,**kwargs):
1091 # TODO: continuation
1092 return ctx.hyper2d({'m+n':[a],'m':[b1],'n':[b2]},
1093 {'m':[c1],'n':[c2]}, x,y, **kwargs)
1095@defun
1096def appellf3(ctx,a1,a2,b1,b2,c,x,y,**kwargs):
1097 outer_polynomial = ctx.isnpint(a1) or ctx.isnpint(b1)
1098 inner_polynomial = ctx.isnpint(a2) or ctx.isnpint(b2)
1099 if not outer_polynomial:
1100 if inner_polynomial or abs(x) > abs(y):
1101 x, y = y, x
1102 a1,a2,b1,b2 = a2,a1,b2,b1
1103 return ctx.hyper2d({'m':[a1,b1],'n':[a2,b2]}, {'m+n':[c]},x,y,**kwargs)
1105@defun
1106def appellf4(ctx,a,b,c1,c2,x,y,**kwargs):
1107 # TODO: continuation
1108 return ctx.hyper2d({'m+n':[a,b]}, {'m':[c1],'n':[c2]},x,y,**kwargs)
1110@defun
1111def hyper2d(ctx, a, b, x, y, **kwargs):
1112 r"""
1113 Sums the generalized 2D hypergeometric series
1115 .. math ::
1117 \sum_{m=0}^{\infty} \sum_{n=0}^{\infty}
1118 \frac{P((a),m,n)}{Q((b),m,n)}
1119 \frac{x^m y^n} {m! n!}
1121 where `(a) = (a_1,\ldots,a_r)`, `(b) = (b_1,\ldots,b_s)` and where
1122 `P` and `Q` are products of rising factorials such as `(a_j)_n` or
1123 `(a_j)_{m+n}`. `P` and `Q` are specified in the form of dicts, with
1124 the `m` and `n` dependence as keys and parameter lists as values.
1125 The supported rising factorials are given in the following table
1126 (note that only a few are supported in `Q`):
1128 +------------+-------------------+--------+
1129 | Key | Rising factorial | `Q` |
1130 +============+===================+========+
1131 | ``'m'`` | `(a_j)_m` | Yes |
1132 +------------+-------------------+--------+
1133 | ``'n'`` | `(a_j)_n` | Yes |
1134 +------------+-------------------+--------+
1135 | ``'m+n'`` | `(a_j)_{m+n}` | Yes |
1136 +------------+-------------------+--------+
1137 | ``'m-n'`` | `(a_j)_{m-n}` | No |
1138 +------------+-------------------+--------+
1139 | ``'n-m'`` | `(a_j)_{n-m}` | No |
1140 +------------+-------------------+--------+
1141 | ``'2m+n'`` | `(a_j)_{2m+n}` | No |
1142 +------------+-------------------+--------+
1143 | ``'2m-n'`` | `(a_j)_{2m-n}` | No |
1144 +------------+-------------------+--------+
1145 | ``'2n-m'`` | `(a_j)_{2n-m}` | No |
1146 +------------+-------------------+--------+
1148 For example, the Appell F1 and F4 functions
1150 .. math ::
1152 F_1 = \sum_{m=0}^{\infty} \sum_{n=0}^{\infty}
1153 \frac{(a)_{m+n} (b)_m (c)_n}{(d)_{m+n}}
1154 \frac{x^m y^n}{m! n!}
1156 F_4 = \sum_{m=0}^{\infty} \sum_{n=0}^{\infty}
1157 \frac{(a)_{m+n} (b)_{m+n}}{(c)_m (d)_{n}}
1158 \frac{x^m y^n}{m! n!}
1160 can be represented respectively as
1162 ``hyper2d({'m+n':[a], 'm':[b], 'n':[c]}, {'m+n':[d]}, x, y)``
1164 ``hyper2d({'m+n':[a,b]}, {'m':[c], 'n':[d]}, x, y)``
1166 More generally, :func:`~mpmath.hyper2d` can evaluate any of the 34 distinct
1167 convergent second-order (generalized Gaussian) hypergeometric
1168 series enumerated by Horn, as well as the Kampe de Feriet
1169 function.
1171 The series is computed by rewriting it so that the inner
1172 series (i.e. the series containing `n` and `y`) has the form of an
1173 ordinary generalized hypergeometric series and thereby can be
1174 evaluated efficiently using :func:`~mpmath.hyper`. If possible,
1175 manually swapping `x` and `y` and the corresponding parameters
1176 can sometimes give better results.
1178 **Examples**
1180 Two separable cases: a product of two geometric series, and a
1181 product of two Gaussian hypergeometric functions::
1183 >>> from mpmath import *
1184 >>> mp.dps = 25; mp.pretty = True
1185 >>> x, y = mpf(0.25), mpf(0.5)
1186 >>> hyper2d({'m':1,'n':1}, {}, x,y)
1187 2.666666666666666666666667
1188 >>> 1/(1-x)/(1-y)
1189 2.666666666666666666666667
1190 >>> hyper2d({'m':[1,2],'n':[3,4]}, {'m':[5],'n':[6]}, x,y)
1191 4.164358531238938319669856
1192 >>> hyp2f1(1,2,5,x)*hyp2f1(3,4,6,y)
1193 4.164358531238938319669856
1195 Some more series that can be done in closed form::
1197 >>> hyper2d({'m':1,'n':1},{'m+n':1},x,y)
1198 2.013417124712514809623881
1199 >>> (exp(x)*x-exp(y)*y)/(x-y)
1200 2.013417124712514809623881
1202 Six of the 34 Horn functions, G1-G3 and H1-H3::
1204 >>> from mpmath import *
1205 >>> mp.dps = 10; mp.pretty = True
1206 >>> x, y = 0.0625, 0.125
1207 >>> a1,a2,b1,b2,c1,c2,d = 1.1,-1.2,-1.3,-1.4,1.5,-1.6,1.7
1208 >>> hyper2d({'m+n':a1,'n-m':b1,'m-n':b2},{},x,y) # G1
1209 1.139090746
1210 >>> nsum(lambda m,n: rf(a1,m+n)*rf(b1,n-m)*rf(b2,m-n)*\
1211 ... x**m*y**n/fac(m)/fac(n), [0,inf], [0,inf])
1212 1.139090746
1213 >>> hyper2d({'m':a1,'n':a2,'n-m':b1,'m-n':b2},{},x,y) # G2
1214 0.9503682696
1215 >>> nsum(lambda m,n: rf(a1,m)*rf(a2,n)*rf(b1,n-m)*rf(b2,m-n)*\
1216 ... x**m*y**n/fac(m)/fac(n), [0,inf], [0,inf])
1217 0.9503682696
1218 >>> hyper2d({'2n-m':a1,'2m-n':a2},{},x,y) # G3
1219 1.029372029
1220 >>> nsum(lambda m,n: rf(a1,2*n-m)*rf(a2,2*m-n)*\
1221 ... x**m*y**n/fac(m)/fac(n), [0,inf], [0,inf])
1222 1.029372029
1223 >>> hyper2d({'m-n':a1,'m+n':b1,'n':c1},{'m':d},x,y) # H1
1224 -1.605331256
1225 >>> nsum(lambda m,n: rf(a1,m-n)*rf(b1,m+n)*rf(c1,n)/rf(d,m)*\
1226 ... x**m*y**n/fac(m)/fac(n), [0,inf], [0,inf])
1227 -1.605331256
1228 >>> hyper2d({'m-n':a1,'m':b1,'n':[c1,c2]},{'m':d},x,y) # H2
1229 -2.35405404
1230 >>> nsum(lambda m,n: rf(a1,m-n)*rf(b1,m)*rf(c1,n)*rf(c2,n)/rf(d,m)*\
1231 ... x**m*y**n/fac(m)/fac(n), [0,inf], [0,inf])
1232 -2.35405404
1233 >>> hyper2d({'2m+n':a1,'n':b1},{'m+n':c1},x,y) # H3
1234 0.974479074
1235 >>> nsum(lambda m,n: rf(a1,2*m+n)*rf(b1,n)/rf(c1,m+n)*\
1236 ... x**m*y**n/fac(m)/fac(n), [0,inf], [0,inf])
1237 0.974479074
1239 **References**
1241 1. [SrivastavaKarlsson]_
1242 2. [Weisstein]_ http://mathworld.wolfram.com/HornFunction.html
1243 3. [Weisstein]_ http://mathworld.wolfram.com/AppellHypergeometricFunction.html
1245 """
1246 x = ctx.convert(x)
1247 y = ctx.convert(y)
1248 def parse(dct, key):
1249 args = dct.pop(key, [])
1250 try:
1251 args = list(args)
1252 except TypeError:
1253 args = [args]
1254 return [ctx.convert(arg) for arg in args]
1255 a_s = dict(a)
1256 b_s = dict(b)
1257 a_m = parse(a, 'm')
1258 a_n = parse(a, 'n')
1259 a_m_add_n = parse(a, 'm+n')
1260 a_m_sub_n = parse(a, 'm-n')
1261 a_n_sub_m = parse(a, 'n-m')
1262 a_2m_add_n = parse(a, '2m+n')
1263 a_2m_sub_n = parse(a, '2m-n')
1264 a_2n_sub_m = parse(a, '2n-m')
1265 b_m = parse(b, 'm')
1266 b_n = parse(b, 'n')
1267 b_m_add_n = parse(b, 'm+n')
1268 if a: raise ValueError("unsupported key: %r" % a.keys()[0])
1269 if b: raise ValueError("unsupported key: %r" % b.keys()[0])
1270 s = 0
1271 outer = ctx.one
1272 m = ctx.mpf(0)
1273 ok_count = 0
1274 prec = ctx.prec
1275 maxterms = kwargs.get('maxterms', 20*prec)
1276 try:
1277 ctx.prec += 10
1278 tol = +ctx.eps
1279 while 1:
1280 inner_sign = 1
1281 outer_sign = 1
1282 inner_a = list(a_n)
1283 inner_b = list(b_n)
1284 outer_a = [a+m for a in a_m]
1285 outer_b = [b+m for b in b_m]
1286 # (a)_{m+n} = (a)_m (a+m)_n
1287 for a in a_m_add_n:
1288 a = a+m
1289 inner_a.append(a)
1290 outer_a.append(a)
1291 # (b)_{m+n} = (b)_m (b+m)_n
1292 for b in b_m_add_n:
1293 b = b+m
1294 inner_b.append(b)
1295 outer_b.append(b)
1296 # (a)_{n-m} = (a-m)_n / (a-m)_m
1297 for a in a_n_sub_m:
1298 inner_a.append(a-m)
1299 outer_b.append(a-m-1)
1300 # (a)_{m-n} = (-1)^(m+n) (1-a-m)_m / (1-a-m)_n
1301 for a in a_m_sub_n:
1302 inner_sign *= (-1)
1303 outer_sign *= (-1)**(m)
1304 inner_b.append(1-a-m)
1305 outer_a.append(-a-m)
1306 # (a)_{2m+n} = (a)_{2m} (a+2m)_n
1307 for a in a_2m_add_n:
1308 inner_a.append(a+2*m)
1309 outer_a.append((a+2*m)*(1+a+2*m))
1310 # (a)_{2m-n} = (-1)^(2m+n) (1-a-2m)_{2m} / (1-a-2m)_n
1311 for a in a_2m_sub_n:
1312 inner_sign *= (-1)
1313 inner_b.append(1-a-2*m)
1314 outer_a.append((a+2*m)*(1+a+2*m))
1315 # (a)_{2n-m} = 4^n ((a-m)/2)_n ((a-m+1)/2)_n / (a-m)_m
1316 for a in a_2n_sub_m:
1317 inner_sign *= 4
1318 inner_a.append(0.5*(a-m))
1319 inner_a.append(0.5*(a-m+1))
1320 outer_b.append(a-m-1)
1321 inner = ctx.hyper(inner_a, inner_b, inner_sign*y,
1322 zeroprec=ctx.prec, **kwargs)
1323 term = outer * inner * outer_sign
1324 if abs(term) < tol:
1325 ok_count += 1
1326 else:
1327 ok_count = 0
1328 if ok_count >= 3 or not outer:
1329 break
1330 s += term
1331 for a in outer_a: outer *= a
1332 for b in outer_b: outer /= b
1333 m += 1
1334 outer = outer * x / m
1335 if m > maxterms:
1336 raise ctx.NoConvergence("maxterms exceeded in hyper2d")
1337 finally:
1338 ctx.prec = prec
1339 return +s
1341"""
1342@defun
1343def kampe_de_feriet(ctx,a,b,c,d,e,f,x,y,**kwargs):
1344 return ctx.hyper2d({'m+n':a,'m':b,'n':c},
1345 {'m+n':d,'m':e,'n':f}, x,y, **kwargs)
1346"""
1348@defun
1349def bihyper(ctx, a_s, b_s, z, **kwargs):
1350 r"""
1351 Evaluates the bilateral hypergeometric series
1353 .. math ::
1355 \,_AH_B(a_1, \ldots, a_k; b_1, \ldots, b_B; z) =
1356 \sum_{n=-\infty}^{\infty}
1357 \frac{(a_1)_n \ldots (a_A)_n}
1358 {(b_1)_n \ldots (b_B)_n} \, z^n
1360 where, for direct convergence, `A = B` and `|z| = 1`, although a
1361 regularized sum exists more generally by considering the
1362 bilateral series as a sum of two ordinary hypergeometric
1363 functions. In order for the series to make sense, none of the
1364 parameters may be integers.
1366 **Examples**
1368 The value of `\,_2H_2` at `z = 1` is given by Dougall's formula::
1370 >>> from mpmath import *
1371 >>> mp.dps = 25; mp.pretty = True
1372 >>> a,b,c,d = 0.5, 1.5, 2.25, 3.25
1373 >>> bihyper([a,b],[c,d],1)
1374 -14.49118026212345786148847
1375 >>> gammaprod([c,d,1-a,1-b,c+d-a-b-1],[c-a,d-a,c-b,d-b])
1376 -14.49118026212345786148847
1378 The regularized function `\,_1H_0` can be expressed as the
1379 sum of one `\,_2F_0` function and one `\,_1F_1` function::
1381 >>> a = mpf(0.25)
1382 >>> z = mpf(0.75)
1383 >>> bihyper([a], [], z)
1384 (0.2454393389657273841385582 + 0.2454393389657273841385582j)
1385 >>> hyper([a,1],[],z) + (hyper([1],[1-a],-1/z)-1)
1386 (0.2454393389657273841385582 + 0.2454393389657273841385582j)
1387 >>> hyper([a,1],[],z) + hyper([1],[2-a],-1/z)/z/(a-1)
1388 (0.2454393389657273841385582 + 0.2454393389657273841385582j)
1390 **References**
1392 1. [Slater]_ (chapter 6: "Bilateral Series", pp. 180-189)
1393 2. [Wikipedia]_ http://en.wikipedia.org/wiki/Bilateral_hypergeometric_series
1395 """
1396 z = ctx.convert(z)
1397 c_s = a_s + b_s
1398 p = len(a_s)
1399 q = len(b_s)
1400 if (p, q) == (0,0) or (p, q) == (1,1):
1401 return ctx.zero * z
1402 neg = (p-q) % 2
1403 def h(*c_s):
1404 a_s = list(c_s[:p])
1405 b_s = list(c_s[p:])
1406 aa_s = [2-b for b in b_s]
1407 bb_s = [2-a for a in a_s]
1408 rp = [(-1)**neg * z] + [1-b for b in b_s] + [1-a for a in a_s]
1409 rc = [-1] + [1]*len(b_s) + [-1]*len(a_s)
1410 T1 = [], [], [], [], a_s + [1], b_s, z
1411 T2 = rp, rc, [], [], aa_s + [1], bb_s, (-1)**neg / z
1412 return T1, T2
1413 return ctx.hypercomb(h, c_s, **kwargs)