def build_hypergeometric_formula(func):
"""
Create a formula object representing the hypergeometric function ``func``.
"""
# We know that no `ap` are negative integers, otherwise "detect poly"
# would have kicked in. However, `ap` could be empty. In this case we can
# use a different basis.
# I'm not aware of a basis that works in all cases.
from sympy import zeros, Matrix, eye
z = Dummy('z')
if func.ap:
afactors = [_x + a for a in func.ap]
bfactors = [_x + b - 1 for b in func.bq]
expr = _x*Mul(*bfactors) - z*Mul(*afactors)
poly = Poly(expr, _x)
n = poly.degree()
basis = []
M = zeros(n)
for k in range(n):
a = func.ap[0] + k
basis += [hyper([a] + list(func.ap[1:]), func.bq, z)]
if k < n - 1:
M[k, k] = -a
M[k, k + 1] = a
B = Matrix(basis)
C = Matrix([[1] + [0]*(n - 1)])
derivs = [eye(n)]
for k in range(n):
derivs.append(M*derivs[k])
l = poly.all_coeffs()
l.reverse()
res = [0]*n
for k, c in enumerate(l):
for r, d in enumerate(C*derivs[k]):
res[r] += c*d
for k, c in enumerate(res):
M[n - 1, k] = -c/derivs[n - 1][0, n - 1]/poly.all_coeffs()[0]
return Formula(func, z, None, [], B, C, M)
else:
# Since there are no `ap`, none of the `bq` can be non-positive
# integers.
basis = []
bq = list(func.bq[:])
for i in range(len(bq)):
basis += [hyper([], bq, z)]
bq[i] += 1
basis += [hyper([], bq, z)]
B = Matrix(basis)
n = len(B)
C = Matrix([[1] + [0]*(n - 1)])
M = zeros(n)
M[0, n - 1] = z/Mul(*func.bq)
for k in range(1, n):
M[k, k - 1] = func.bq[k - 1]
M[k, k] = -func.bq[k - 1]
return Formula(func, z, None, [], B, C, M)
评论列表
文章目录
