problems with hypergeometric and maxima_calculus.hgfpoly

[Dima Pasechnik]

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)

---------------------------------------------------------------------------

ValueError                                Traceback (most recent call last)

<ipython-input-4-0bef797c6f48> in <module>()

----> 1 hypergeometric([-Integer(2),-Integer(1)],[Integer(2)],-Integer(1)).n(Integer(100))

/home/scratch/dimpase/sage/sage6.3/local/lib/python2.7/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression._numerical_approx (build/cythonized/sage/symbolic/expression.cpp:26972)()

/home/scratch/dimpase/sage/sage6.3/local/lib/python2.7/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression._convert (build/cythonized/sage/symbolic/expression.cpp:7667)()

/home/scratch/dimpase/sage/sage6.3/local/lib/python2.7/site-packages/sage/functions/hypergeometric.pyc in _evalf_(self, a, b, z, parent, algorithm)

    294         aa = [rational_param_as_tuple(c) for c in a]

    295         bb = [rational_param_as_tuple(c) for c in b]

--> 296         return mpmath_utils.call(hyper, aa, bb, z, parent=parent)

    297 

    298     def _tderivative_(self, a, b, z, *args, **kwargs):

/home/scratch/dimpase/sage/sage6.3/local/lib/python2.7/site-packages/sage/libs/mpmath/utils.so in sage.libs.mpmath.utils.call (build/cythonized/sage/libs/mpmath/utils.c:6309)()

/home/scratch/dimpase/sage/sage6.3/local/lib/python2.7/site-packages/mpmath/functions/hypergeometric.pyc in hyper(ctx, a_s, b_s, z, **kwargs)

    223         elif q == 0: return ctx._hyp1f0(a_s[0][0], z)

    224     elif p == 2:

--> 225         if   q == 1: return ctx._hyp2f1(a_s, b_s, z, **kwargs)

    226         elif q == 2: return ctx._hyp2f2(a_s, b_s, z, **kwargs)

    227         elif q == 3: return ctx._hyp2f3(a_s, b_s, z, **kwargs)

/home/scratch/dimpase/sage/sage6.3/local/lib/python2.7/site-packages/mpmath/functions/hypergeometric.pyc in _hyp2f1(ctx, a_s, b_s, z, **kwargs)

    442     if absz <= 0.8 or (ctx.isint(a) and a <= 0 and a >= -1000) or \

    443                       (ctx.isint(b) and b <= 0 and b >= -1000):

--> 444         return ctx.hypsum(2, 1, (atype, btype, ctype), [a, b, c], z, **kwargs)

    445 

    446     orig = ctx.prec

/home/scratch/dimpase/sage/sage6.3/local/lib/python2.7/site-packages/mpmath/ctx_mp.pyc in hypsum(ctx, p, q, flags, coeffs, z, accurate_small, **kwargs)

    686         while 1:

    687             if extraprec > maxprec:

--> 688                 raise ValueError(ctx._hypsum_msg % (prec, prec+extraprec))

    689             wp = prec + extraprec

    690             if magnitude_check:

ValueError: hypsum() failed to converge to the requested 100 bits of accuracy

using a working precision of 7135 bits. Try with a higher maxprec,

maxterms, or set zeroprec.

sage: 

(and if I try .n(10000) the error is the same, only the last bit reads as

ValueError: hypsum() failed to converge to the requested 10000 bits of accuracy

using a working precision of 66315 bits. Try with a higher maxprec,

maxterms, or set zeroprec.)

and

sage: maxima_calculus.hgfpoly([1,2],[2],-1)

---------------------------------------------------------------------------

TypeError                                 Traceback (most recent call last)

<ipython-input-6-a307d86456bf> in <module>()

----> 1 maxima_calculus.hgfpoly([Integer(1),Integer(2)],[Integer(2)],-Integer(1))

/home/scratch/dimpase/sage/sage6.3/local/lib/python2.7/site-packages/sage/interfaces/interface.pyc in __call__(self, *args, **kwds)

    561 

    562     def __call__(self, *args, **kwds):

--> 563         return self._parent.function_call(self._name, list(args), kwds)

    564 

    565     def _sage_doc_(self):

/home/scratch/dimpase/sage/sage6.3/local/lib/python2.7/site-packages/sage/interfaces/interface.pyc in function_call(self, function, args, kwds)

    487                                        [s.name() for s in args],

    488                                        ['%s=%s'%(key,value.name()) for key, value in kwds.items()])

--> 489         return self.new(s)

    490 

    491     def _function_call_string(self, function, args, kwds):

/home/scratch/dimpase/sage/sage6.3/local/lib/python2.7/site-packages/sage/interfaces/interface.pyc in new(self, code)

    262 

    263     def new(self, code):

--> 264         return self(code)

    265 

    266     ###################################################################

/home/scratch/dimpase/sage/sage6.3/local/lib/python2.7/site-packages/sage/interfaces/interface.pyc in __call__(self, x, name)

    197 

    198         if isinstance(x, basestring):

--> 199             return cls(self, x, name=name)

    200         try:

    201             return self._coerce_from_special_method(x)

/home/scratch/dimpase/sage/sage6.3/local/lib/python2.7/site-packages/sage/interfaces/interface.pyc in __init__(self, parent, value, is_name, name)

    624                 self._name = parent._create(value, name=name)

    625             except (TypeError, RuntimeError, ValueError) as x:

--> 626                 raise TypeError(x)

    627 

    628     def _latex_(self):

TypeError: ECL says: In function -, the value of the only argument is

  NIL

which is not of the expected type NUMBER

sage: 

-----------------------------------------------------------------------------

is any of these known? In the 1st case, looks like numeric evaluation of qFp by mpmath is

seriously broken...

[Fredrik Johansson]

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

[kcrisman]

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.

Except then Sage would have to know when it is a polynomial, and probably we would need to ask Maxima for that (assuming it knows).  So maybe not completely trivial to make sure it works.

- kcrisman 

[Fredrik Johansson]

It's a polynomial when any of the first parameters is a nonpositive integer.

Fredrik

[kcrisman]

>> 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.
>>
>
> Except then Sage would have to know when it is a polynomial, and probably we
> would need to ask Maxima for that (assuming it knows).  So maybe not
> completely trivial to make sure it works.

It's a polynomial when any of the first parameters is a nonpositive integer.

Is that "if and only if"?  That would certainly be convenient. 

[John Cremona]

I believe this is true for the same reason as for the power series
expansion of (1-x)^{-n}, there are rising factorials and these hit
zero if and only if they start at a nonpositive integer. But you
should check the definition of hypergeometric function to be 100%
sure.

John

>
> --
> 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 sage-devel+...@googlegroups.com.
> To post to this group, send email to sage-...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sage-devel.
> For more options, visit https://groups.google.com/d/optout.

[Dima Pasechnik]

isn't it more or less directly by definition ?

Dima

[Fredrik Johansson]

Yes, that follows immediately from the definition.

I should clarify that the hypergeometric function can be zero for
rational input even if it is not a polynomial in z. For example
1F1(3/2, 1/2, z) = (1 + 2*z) * exp(z) (you can verify that mpmath does
not converge with z = -0.5). You are less likely to run into such
zeros, though.

Fredrik

[kcrisman]

>> >> 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.
>> >>
>> >
>> > Except then Sage would have to know when it is a polynomial, and
>> > probably we
>> > would need to ask Maxima for that (assuming it knows).  So maybe not
>> > completely trivial to make sure it works.
>>
>> It's a polynomial when any of the first parameters is a nonpositive
>> integer.
>>
>
> Is that "if and only if"?  That would certainly be convenient.

Yes, that follows immediately from the definition.

Well, then I guess this shouldn't be too hard to implement... I will admit that I will not be the one doing it, but sounds like a no-brainer. 

[Dima Pasechnik]

IMHO it would be great to have this in mpmath rather than a Sage workaround.
(AFAIK, hypergeometric_U is vulnerable to the very same problem, so this
would mean the workaround would either be needed in several places...)

Dima

[Nils Bruin]

On Friday, September 26, 2014 6:55:39 AM UTC-7, Dima Pasechnik wrote:

IMHO it would be great to have this in mpmath rather than a Sage workaround.
(AFAIK, hypergeometric_U is vulnerable to the very same problem, so this
would mean the workaround would either be needed in several places...)

It would be good if mpmath used a different approach when it's known a hypergeometric function is a polynomial, but it's a little orthogonal to what happens in sage. By the time we hit mpmath, we're doing multiprecision float computation. Polynomial hypergeometric functions can easily be evaluated at arbitrary inputs (rational, algebraic, symbolic). We should not be bothering with multiprecision floats if we don't have to, or at least only do so if there's a benefit and we can control precision issues.

[Ralf Stephan]

On Friday, September 26, 2014 4:25:58 PM UTC+2, Nils Bruin wrote:

It would be good if mpmath used a different approach when it's known a hypergeometric function is a polynomial, but it's a little orthogonal to what happens in sage. By the time we hit mpmath, we're doing multiprecision float computation. Polynomial hypergeometric functions can easily be evaluated at arbitrary inputs (rational, algebraic, symbolic). We should not be bothering with multiprecision floats if we don't have to, or at least only do so if there's a benefit and we can control precision issues.

In most symbolic special functions we check for special values and, if found, evaluate immediately. This is not so in hypergeometric.py but would be trivial to add. I can only speculate that the reason was the then forced call of maxima. Since this doesn't seem a problem with commenters please review

http://trac.sagemath.org/ticket/17066

and this can be simply emulated until inclusion of the ticket by saying

sage: hypergeometric([-2,-1],[2],-1).simplify_hypergeometric()
0

No need to explicitly call Maxima here.

Regards,