Ask Your Question
1

Hypergeometric Bug?

asked 5 years ago

rrogers gravatar image

updated 5 years ago

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.

Preview: (hide)

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 ( 5 years ago )

2 Answers

Sort by » oldest newest most voted
1

answered 5 years ago

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.

Preview: (hide)
link

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 ( 5 years ago )
1

answered 5 years ago

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.

Preview: (hide)
link

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: 5 years ago

Seen: 439 times

Last updated: Oct 21 '19