# Hypergeometric Bug?

The following code shows a (typical ?) failure. F([m,a,b;m,c;z) should and does for positive m; cancel m. but when m is negative "0" is declared prematurely

--Juptyter code

def pochhammer(x,a):
return rising_factorial(x,a)

def fn2ref(a,b,c,z,m): return hypergeometric([a,b],[c],z)

def fn2lm(a,b,c,z,m):
return hypergeometric([m,a,b],[c,m],z)

def fn2l(a,b,c,z,m):
return hypergeometric([-m,a,b],[c,-m],z)

def fn1l(a,b,c,z,m):
var('k')
return sum((pochhammer(a,k)*pochhammer(b,k)/pochhammer(c,k))*(z^k/factorial(k)),k,0,m)

print "2F1(a,b;c;z,5) = ", fn2ref(1,2,3,.5,5)
print "3F2(m,a,b;m,c;z) = ", fn2lm(1,2,3,.5,5)
print "3F2(-m,a,b;-m,c;z) = ", fn2l(1,2,3,.5,5)
print "m terms of 2F1(a,b;c;z) = ", fn1l(1,2,3,.5,5)


-- produces
2F1(a,b;c;z,5) = 1.54517744447956
3F2(m,a,b;m,c;z) = 1.54517744447956
3F2(-m,a,b;-m,c;z) = 1.53809523809524
m terms of 2F1(a,b;c;z) = 1.538095238095238

--Whereas Maxima code

fn2ref(a,b,c,z):=hypergeometric([a,b],[c],z);
fn2lm(a,b,c,z,m)      :=hypergeometric([a,b,m],[c,m],z);
fn2l(a,b,c,z,m)      :=hypergeometric([a,b,-m],[c,-m],z);
fn1l(a,b,c,z,m):=
sum((pochhammer(a,k)*pochhammer(b,k)/(pochhammer(c,k))*z^k/factorial(k)),k,0,m);

print ("2F1(a,b;c;z) = ",fn2ref(1,2,3,.5))$print("3F2(m,a,b;m,c;z) = ",fn2lm(1,2,3,.5,5))$

print("3F2(-m,a,b;-m,c;z) = ",fn2l(1,2,3,.5,5))$print("m terms of 2F1(a,b;c;z) = ",fn1l(1,2,3,.5,5))$


-- Produces
2F1(a,b;c;z) = 1.545177444479563
3F2(m,a,b;m,c;z) = 1.545177444479563
3F2(-m,a,b;-m,c;z) = 1.545177444479563
m terms of 2F1(a,b;c;z) = 1.538095238095238

--
Which accords with the definition and Luke "The Special Functions and their Approximations" which is a pretty standard.

edit retag close merge delete

To display blocks of code or error messages, skip a line above and below, and do one of the following (all give the same result):

• indent all code lines with 4 spaces
• select all code lines and click the "code" button (the icon with '101 010')
• select all code lines and hit ctrl-K

For instance, typing

If we define f by

def f(x, y, z):
return x * y * z

then f(2, 3, 5) returns 30 but f(2*3*5) gives:

TypeError: f() takes exactly 3 arguments (1 given)


produces:

If we define f by

def f(x, y, z):
return x * y * z


then f(2, 3, 5) returns 30 but f(2*3*5) gives:

TypeError: f() takes exactly 3 arguments (1 given)


( 2019-10-06 10:32:16 -0500 )edit

Sort by » oldest newest most voted

The reason given in the documentation cited in the comment above is debatable. Read this section in DLMF https://dlmf.nist.gov/15.2#ii which deals with a similar case and says: "Because of the analytic properties [...] it is usually legitimate to take limits in formulas involving functions that are undefined for certain values of the parameters."

For instance Maple gives the same values as Maxima.

more

Sage uses mpmath for numerical evaluation. The documentation for the underlying function says explicitly:

If a nonpositive integer −n appears in both a_s and b_s, this parameter cannot be unambiguously removed since it creates a term 0 / 0. In this case the hypergeometric series is understood to terminate before the division by zero occurs. This convention is consistent with Mathematica.

I have confirmed that Mathematica gives the same answers as Sage. Apparently Maxima chooses a different convention, so be aware of this.

more

Really?? (sigh) this type of thing is why I always double-check symbolic/numeric algebra packages (:(: and have to write "shells" (incidently it's wrong).

( 2019-10-12 12:21:10 -0500 )edit