always simplify hypergeometric() when it's a polynomial

[rws]

See ​https://groups.google.com/forum/?hl=en#!topic/sage-devel/8z2HyH5Y8Dc with Fredrik's assessment:

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

[rws]

This patch just changes hypergeometric() to always invoke Maxima_lib.hgfred().

---

New commits:

​b6f0fbc

17066: always try to simplify

[kcrisman]

Maybe document the problem at ask.sagemath is fixed, too?

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

was the example Dima posted on sage-devel...

[rws]

The "solution" is tongue-in-cheek, anyway, at least its application should be restricted in order to not burn CPU needlessly with higher orders of the function.

[kcrisman]

Perhaps, but it wouldn't be hard to restrict it properly, right?

[kcrisman]

This is related to the following ​ask.sagemath question as well.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​661487a	Merge branch 'develop' into t/17066/work_around_mpmath_problem_with_hypergeometric___zeroes
​dff348c	17066: fix/use closed_form() instead of Maxima

[rws]

Replying to kcrisman:

This is related to the following ​ask.sagemath question as well.

Indeed. While what you're saying there about decidability is certainly right, we can bit off a good chunk using known closed forms. I have now fixed a bug in closed_form() and it looks like using that Sage function (provided by Fredrik with the hypergeometric code) instead of Maxima is not only faster but gives also better results:

sage: hypergeometric([3/2],[1/2,-1/2],5)        # using closed_form()
-2*sqrt(5)*sinh(2*sqrt(5)) - 19*cosh(2*sqrt(5))
sage: from sage.interfaces.maxima_lib import maxima_lib, max_to_sr
sage: max_to_sr(maxima_lib.hgfred([3/2],[1/2,-1/2],5).ecl())
-4*(5*5^(1/4))*sqrt(1/5)*sqrt(pi)*sqrt(sqrt(5)/pi)*cosh(2*sqrt(5)) - 2*5^(3/4)*sqrt(pi)*bessel_I(-3/2, 2*sqrt(5))

The latter, while numerically correct, is difficult to simplify further.

There are several doctests failing where I need to decide how to fix or rewrite them.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​a225383

17066: restrict to nonpositive first arg; fix doctests

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​66ed0c0

17066: more doctests

[rws]

Had to restrict the simplification, after all, as advised.

[kcrisman]

Yeah, make sure nothing is wrong!

I know that I won't be able to finish the review of this immediately. But I do notice that all the hold=True will be quite confusing in the doctests. I recommend that we add some more doctests that test the methods in question with non-simplifying stuff, and then make it clear to the reader why we are all of a sudden using hold=True. On the plus side, hold=True seems to have worked great!

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​55daf81

17066: explain why we add hold=True to some doctests

[rws]

Does this take #10169 comment 4 in account?

[mmezzarobba]

Expanding hypergeometric functions in polynomial closed form before evaluating them can lead to numerical cancellation. Compare for example

sage: hypergeometric([1,-100],[3], 0.7)
0.0278923450568346
sage: hypergeometric([1,-100],[3], Reals(1000)(0.7)).n()
0.0278923450568346
sage: sage.functions.hypergeometric.closed_form(hypergeometric([1,-100],[3],x)).subs(x=0.7)
-376.859920391941

Won't that be a problem with this ticket? (I didn't read beyond the ticket's description, so please ignore my comment if it is beside the point.)

[rws]

Replying to mmezzarobba:

Won't that be a problem with this ticket?

It appears the only solutions are (apart from relying on the user to supply needed precision): A) restrict the range of parameters where ticket is applied B) do not simplify; just fix the evaluation problem C) declare the problem to be in n()

I mean really, n() makes bold claims:

   Returns a numerical approximation of an object "x" with at least
   "prec" bits (or decimal "digits") of precision.

but won't deliver when x is a polynomial?

[fredrik.johansson]

When it's a polynomial of low degree and the inputs are not too large, you could do everything with rational arithmetic and convert to a correctly rounded floating-point number in the end. Or you could use interval arithmetic.

[mmezzarobba]

Replying to rws:

It appears the only solutions are A) restrict the range of parameters where ticket is applied B) do not simplify; just fix the evaluation problem C) rely on the user to supply needed precision; add warning

If (as you seem to say in your post to sage-devel) the result of the ticket is that

sage: hypergeometric([1,-100],[3],x)

returns a polynomial while

sage: hypergeometric([1,-100],[3],0.7)

still calls a numerically stable evaluation algorithm, then I don't think there is a problem. (Ideally, hypergeometric([1,-100],[3],x,hold=True).subs(x=0.7) should also call the numerically stable evaluation algorithm, though.)

What I was concerned about was that the implementation of n() for an "atomic-looking" expression with dedicated numerical evaluation code would end up using a numerically unstable formula internally without adapting the precision or warning the user. I'm not saying such a behaviour is an absolute no-go, but I do think it is something we should avoid when possible.

[rws]

It's not that I don't see your point. The real fix for this cannot be provided by hypergeometric though but must be done in subs, if at all.

[fredrik.johansson]

Since mpmath uses a numerically stable algorithm but fails on zero, you could potentially also try calling mpmath, catch the exception, and then check symbolically if it's exactly zero.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​290b47f	Merge branch 'develop' into t/17066/work_around_mpmath_problem_with_hypergeometric___zeroes
​0e0fef6	17066: add warning

[kcrisman]

Expanding hypergeometric functions in polynomial closed form before evaluating them can lead to numerical cancellation. Compare for example

sage: hypergeometric([1,-100],[3], 0.7)
0.0278923450568346
sage: hypergeometric([1,-100],[3], Reals(1000)(0.7)).n()
0.0278923450568346
sage: sage.functions.hypergeometric.closed_form(hypergeometric([1,-100],[3],x)).subs(x=0.7)
-376.859920391941

Won't that be a problem with this ticket? (I didn't read beyond the ticket's description, so please ignore my comment if it is beside the point.)

Just so I follow, what you are saying is that the polynomial (rational function, I guess) given will be correct, but that substituting floating point numbers are likely to be inaccurate unless carried out to a very high degree of precision? I'm torn here because on the one hand, presumably that could happen with any rational function, there are standard examples given in a lot of calculus books, but on the other hand we "expect" hg to return a hg function and not a polynomial.

The real fix for this cannot be provided by hypergeometric though but must be done in subs, if at all.

Well, it could provide a good reason not to "always simplify hypergeometric when it's a polynomial"... I didn't realize that people meant rational function when they said polynomial. Is it possible to really only simplify when it's a true polynomial, or is that a very unlikely event?

[mmezzarobba]

Replying to kcrisman:

Just so I follow, what you are saying is that the polynomial (rational function, I guess)

No, it is really a polynomial—but one of degree 100, with large coefficients that alternate in sign.

given will be correct, but that substituting floating point numbers are likely to be inaccurate unless carried out to a very high degree of precision?

Precisely. (Well, not necessarily a very high precision, but probably something like 80-100 guard bits in my example.)

[kcrisman]

Just so I follow, what you are saying is that the polynomial (rational function, I guess)

No, it is really a polynomial—but one of degree 100, with large coefficients that alternate in sign.

Oh, I misread some of the output in the example I tested - hard to see where the dividing signs go :-) So it's the alternation that is the cancellation problem, got it.

given will be correct, but that substituting floating point numbers are likely to be inaccurate unless carried out to a very high degree of precision?

Precisely. (Well, not necessarily a very high precision, but probably something like 80-100 guard bits in my example.)

Anyway, 53 default yields madness, as we see in comment:16. Not good, but I'm not sure what the "right" answer is; I'm leaning toward saying we should not do this automatically, then, but making the simplification a lot easier to find (which is probably already done in this ticket). There are likely as many people using hg functions for numerics as for symbolics, right?

[rws]

sage -t --long src/sage/symbolic/expression.pyx  # 4 doctests failed
sage -t --long src/sage/functions/bessel.py  # 1 doctest failed
sage -t --long src/sage/functions/hypergeometric.py  # 102 doctests failed