mpmath / mpmath

Python library for arbitrary-precision floating-point arithmetic
http://mpmath.org
BSD 3-Clause "New" or "Revised" License
966 stars 184 forks source link

meijerg(((0.5, 1, 1), ()), ((0.5, 1), (0,)), 1.0, 1) returns nan #378

Open skirpichev opened 6 years ago

skirpichev commented 6 years ago

Unless I miss something, this should be MeijerG[{{1/2, 1, 1}, {}}, {{1/2, 1}, {0}}, 1] ~ 4.355 in the Mathematica.

fredrik-johansson commented 6 years ago

I did not check this thoroughly, but it looks like a special case where the Meijer G-function is a limiting sum of two hypergeometric functions that both blow up at z = 1. The hypercomb code handles limits with respect to parameters automatically, but not with respect to the argument. I don't know if this can be fixed in a sensible way in hypercomb or if a special case for the Meijer G-function is required.

There is a simple enough workaround when you know that this case is occurring:

>>> with extraprec(mp.prec): y = meijerg(((0.5, 1, 1), ()), ((0.5, 1), (0,)), 1 - sqrt(eps))
... 
>>> y
mpf('4.355172180607204')
rouckas commented 5 years ago

This is probably related: mp.meijerg(((), (1/2,)), ((-3/2, 0), ()), 1/2) fails with

ValueError: hypsum() failed to converge to the requested 63 bits of accuracy
using a working precision of 3578 bits. Try with a higher maxprec,
maxterms, or set zeroprec.

adding epsilon to the argument also works as a workaround.

skirpichev commented 2 months ago

@rouckas, did you try suggestion "Try with a higher maxprec, maxterms, or set zeroprec"?

rouckas commented 2 months ago

@skirpichev now I tried, and indeed, setting zeroprec works:

In [1]: mp.meijerg(((), (1/2,)), ((-3/2, 0), ()), 1/2, zeroprec=mp.prec)
Out[1]: mpf('1.2130613194252668')

I think I didn't try these options before because they weren't documented. But looking into the source code, it seems zeroprec basically determines the absolute precision (relative to unity), so it makes sense now.