Ask Your Question
1

Hypergeometric Bug?

asked 2019-10-06 16:16:19 +0100

rrogers gravatar image

updated 2019-10-14 13:45:17 +0100

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 flag offensive close merge delete

Comments

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)

Please edit your question to do that.

slelievre gravatar imageslelievre ( 2019-10-06 17:32:16 +0100 )edit

2 Answers

Sort by » oldest newest most voted
1

answered 2019-10-08 22:57:18 +0100

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.

edit flag offensive delete link more

Comments

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).

rrogers gravatar imagerrogers ( 2019-10-12 19:21:10 +0100 )edit
1

answered 2019-10-21 00:07:22 +0100

Peter Luschny gravatar image

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.

edit flag offensive delete link more

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2019-10-06 16:16:19 +0100

Seen: 424 times

Last updated: Oct 21 '19