On Wednesday, September 24, 2014 12:35:50 PM UTC+2, Dima Pasechnik wrote:
>
> there is quite a bit of weirdness mentioned here:
>
> http://ask.sagemath.org/question/24257/a-hypergeometric-series/
>
> e.g.
>
> sage: hypergeometric([-2,-1],[2],-1).n(100)
> ---------------------------------------------------------------------------
>
>
> -----------------------------------------------------------------------------
>
> is any of these known? In the 1st case, looks like numeric evaluation of
> qFp by mpmath is
> seriously broken...
>
>
hyp2f1(-2,-1,2,-1) doesn't converge because the value is zero, and the
hypergeometric function evaluation code in mpmath always tries to achieve
full *relative* accuracy. One solution is to do:
>>> hyp2f1(-2,-1,2,-1,zeroprec=1000)
mpf('0.0')
This is saying that if the result cannot be distinguished from zero when
using 1000 bits of working precision, it is probably actually zero and not
just merely something very small (so don't be fussy about it).
In fact, hyp2f1(-2,-1,2,z) is just the polynomial 1+z. mpmath has no way to
distinguish 1 + (-1) from 1 + (-1 + eps) in which the eps disappears due
rounding, so it has to assume that something may have disappeared. Indeed:
>>> mp.dps = 1000
>>> z = mpf(-1) + mpf("1e-900")
>>> mp.dps = 15
>>> hyper([-2,-1],[2],z) # conservatively throws an error
Traceback (most recent call last):
...
ValueError: hypsum() failed to converge to the requested 53 bits of accuracy
using a working precision of 3568 bits. Try with a higher maxprec,
maxterms, or set zeroprec.
>>> hyper([-2,-1],[2],z,zeroprec=1000) # WRONG!!! (relatively speaking)
mpf('0.0')
>>> hyper([-2,-1],[2],z,maxprec=10000) # correct
mpf('9.9999999999999998e-901')
Exact zero detection could be done at least in trivial cases by identifying
easy factors such as z-1 or z+1. In general, I think one basically has to
fall back to using exact rational arithmetic for the evaluation. This would
be worth implementing in mpmath, because the current behavior is definitely
annoying.
For Sage, fixing the problem is actually trivial: when the hypergeometric
function is a polynomial (and at least when the inputs are exact), don't
call mpmath; just evaluate the polynomial directly and then call .n() on
the result.
Fredrik
--
You received this message because you are subscribed to the Google Groups
"sage-devel" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To post to this group, send email to [email protected].
Visit this group at http://groups.google.com/group/sage-devel.
For more options, visit https://groups.google.com/d/optout.
