Sage: Ticket #17066: always simplify hypergeometric() when it's a polynomial
https://trac.sagemath.org/ticket/17066
<p>
See <a class="ext-link" href="https://groups.google.com/forum/?hl=en#!topic/sage-devel/8z2HyH5Y8Dc"><span class="icon"></span>https://groups.google.com/forum/?hl=en#!topic/sage-devel/8z2HyH5Y8Dc</a> with Fredrik's assessment:
</p>
<p>
"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."
</p>
en-usSagehttps://trac.sagemath.org/chrome/site/logo_sagemath_trac.png
https://trac.sagemath.org/ticket/17066
Trac 1.1.6rwsTue, 30 Sep 2014 16:47:12 GMTbranch set
https://trac.sagemath.org/ticket/17066#comment:1
https://trac.sagemath.org/ticket/17066#comment:1
<ul>
<li><strong>branch</strong>
set to <em>u/rws/work_around_mpmath_problem_with_hypergeometric___zeroes</em>
</li>
</ul>
TicketrwsTue, 30 Sep 2014 16:50:14 GMTstatus, summary changed; commit, author set
https://trac.sagemath.org/ticket/17066#comment:2
https://trac.sagemath.org/ticket/17066#comment:2
<ul>
<li><strong>status</strong>
changed from <em>new</em> to <em>needs_review</em>
</li>
<li><strong>commit</strong>
set to <em>b6f0fbc3983b520c28330e49ec38d79aadf1796f</em>
</li>
<li><strong>summary</strong>
changed from <em>work around mpmath problem with hypergeometric() zeroes</em> to <em>always try to simplify hypergeometric()</em>
</li>
<li><strong>author</strong>
set to <em>Ralf Stephan</em>
</li>
</ul>
<p>
This patch just changes <code>hypergeometric()</code> to always invoke <code>Maxima_lib.hgfred()</code>.
</p>
<hr />
<p>
New commits:
</p>
<table class="wiki">
<tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=b6f0fbc3983b520c28330e49ec38d79aadf1796f"><span class="icon"></span>b6f0fbc</a></td><td><code>17066: always try to simplify</code>
</td></tr></table>
TicketkcrismanWed, 01 Oct 2014 18:10:52 GMTstatus changed
https://trac.sagemath.org/ticket/17066#comment:3
https://trac.sagemath.org/ticket/17066#comment:3
<ul>
<li><strong>status</strong>
changed from <em>needs_review</em> to <em>needs_work</em>
</li>
</ul>
<p>
Maybe document the problem at ask.sagemath is fixed, too?
</p>
<pre class="wiki">hypergeometric([-2,-1],[2],-1).n(100)
</pre><p>
was the example Dima posted on sage-devel...
</p>
TicketrwsThu, 02 Oct 2014 09:41:00 GMT
https://trac.sagemath.org/ticket/17066#comment:4
https://trac.sagemath.org/ticket/17066#comment:4
<p>
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.
</p>
TicketkcrismanThu, 02 Oct 2014 10:16:41 GMT
https://trac.sagemath.org/ticket/17066#comment:5
https://trac.sagemath.org/ticket/17066#comment:5
<p>
Perhaps, but it wouldn't be hard to restrict it properly, right?
</p>
TicketkcrismanWed, 29 Oct 2014 16:09:44 GMT
https://trac.sagemath.org/ticket/17066#comment:6
https://trac.sagemath.org/ticket/17066#comment:6
<p>
This is related to the following <a class="ext-link" href="http://ask.sagemath.org/question/24677/a-simple-hypergeometric-function-fails/"><span class="icon"></span>ask.sagemath question</a> as well.
</p>
TicketgitFri, 31 Oct 2014 09:23:14 GMTcommit changed
https://trac.sagemath.org/ticket/17066#comment:7
https://trac.sagemath.org/ticket/17066#comment:7
<ul>
<li><strong>commit</strong>
changed from <em>b6f0fbc3983b520c28330e49ec38d79aadf1796f</em> to <em>dff348c1c81579de07e88485e05d5e40b0d31942</em>
</li>
</ul>
<p>
Branch pushed to git repo; I updated commit sha1. New commits:
</p>
<table class="wiki">
<tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=661487a058af4d581ff87b43d7fed55ecb494933"><span class="icon"></span>661487a</a></td><td><code>Merge branch 'develop' into t/17066/work_around_mpmath_problem_with_hypergeometric___zeroes</code>
</td></tr><tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=dff348c1c81579de07e88485e05d5e40b0d31942"><span class="icon"></span>dff348c</a></td><td><code>17066: fix/use closed_form() instead of Maxima</code>
</td></tr></table>
TicketrwsFri, 31 Oct 2014 09:23:56 GMT
https://trac.sagemath.org/ticket/17066#comment:8
https://trac.sagemath.org/ticket/17066#comment:8
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/17066#comment:6" title="Comment 6">kcrisman</a>:
</p>
<blockquote class="citation">
<p>
This is related to the following <a class="ext-link" href="http://ask.sagemath.org/question/24677/a-simple-hypergeometric-function-fails/"><span class="icon"></span>ask.sagemath question</a> as well.
</p>
</blockquote>
<p>
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 <code>closed_form()</code> 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:
</p>
<pre class="wiki">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))
</pre><p>
The latter, while numerically correct, is difficult to simplify further.
</p>
<p>
There are several doctests failing where I need to decide how to fix or rewrite them.
</p>
TicketgitFri, 31 Oct 2014 10:21:55 GMTcommit changed
https://trac.sagemath.org/ticket/17066#comment:9
https://trac.sagemath.org/ticket/17066#comment:9
<ul>
<li><strong>commit</strong>
changed from <em>dff348c1c81579de07e88485e05d5e40b0d31942</em> to <em>a2253833d3edd9cd6fa9ccfc9509592ad2b9fad4</em>
</li>
</ul>
<p>
Branch pushed to git repo; I updated commit sha1. New commits:
</p>
<table class="wiki">
<tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=a2253833d3edd9cd6fa9ccfc9509592ad2b9fad4"><span class="icon"></span>a225383</a></td><td><code>17066: restrict to nonpositive first arg; fix doctests</code>
</td></tr></table>
TicketgitFri, 31 Oct 2014 10:31:40 GMTcommit changed
https://trac.sagemath.org/ticket/17066#comment:10
https://trac.sagemath.org/ticket/17066#comment:10
<ul>
<li><strong>commit</strong>
changed from <em>a2253833d3edd9cd6fa9ccfc9509592ad2b9fad4</em> to <em>66ed0c0586f2ec1b78a28e76a1f34f258b6f4f96</em>
</li>
</ul>
<p>
Branch pushed to git repo; I updated commit sha1. New commits:
</p>
<table class="wiki">
<tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=66ed0c0586f2ec1b78a28e76a1f34f258b6f4f96"><span class="icon"></span>66ed0c0</a></td><td><code>17066: more doctests</code>
</td></tr></table>
TicketrwsFri, 31 Oct 2014 10:32:49 GMTstatus changed
https://trac.sagemath.org/ticket/17066#comment:11
https://trac.sagemath.org/ticket/17066#comment:11
<ul>
<li><strong>status</strong>
changed from <em>needs_work</em> to <em>needs_review</em>
</li>
</ul>
<p>
Had to restrict the simplification, after all, as advised.
</p>
TicketkcrismanFri, 31 Oct 2014 12:10:30 GMT
https://trac.sagemath.org/ticket/17066#comment:12
https://trac.sagemath.org/ticket/17066#comment:12
<p>
Yeah, make sure nothing is <em>wrong</em>!
</p>
<p>
I know that I won't be able to finish the review of this immediately. But I do notice that all the <code>hold=True</code> 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 <code>hold=True</code>. On the plus side, <code>hold=True</code> seems to have worked great!
</p>
TicketgitSat, 01 Nov 2014 15:48:41 GMTcommit changed
https://trac.sagemath.org/ticket/17066#comment:13
https://trac.sagemath.org/ticket/17066#comment:13
<ul>
<li><strong>commit</strong>
changed from <em>66ed0c0586f2ec1b78a28e76a1f34f258b6f4f96</em> to <em>55daf81080f46c54429e629b65e8049acfdb6d00</em>
</li>
</ul>
<p>
Branch pushed to git repo; I updated commit sha1. New commits:
</p>
<table class="wiki">
<tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=55daf81080f46c54429e629b65e8049acfdb6d00"><span class="icon"></span>55daf81</a></td><td><code>17066: explain why we add hold=True to some doctests</code>
</td></tr></table>
TicketrwsSat, 01 Nov 2014 15:50:06 GMTsummary changed
https://trac.sagemath.org/ticket/17066#comment:14
https://trac.sagemath.org/ticket/17066#comment:14
<ul>
<li><strong>summary</strong>
changed from <em>always try to simplify hypergeometric()</em> to <em>always simplify hypergeometric() when it's a polynomial</em>
</li>
</ul>
TicketrwsThu, 25 Dec 2014 09:07:02 GMT
https://trac.sagemath.org/ticket/17066#comment:15
https://trac.sagemath.org/ticket/17066#comment:15
<p>
Does this take <a class="new ticket" href="https://trac.sagemath.org/ticket/10169" title="defect: Operands and Operator of symbolic expressions (new)">#10169</a> comment 4 in account?
</p>
TicketmmezzarobbaSun, 18 Jan 2015 19:34:07 GMT
https://trac.sagemath.org/ticket/17066#comment:16
https://trac.sagemath.org/ticket/17066#comment:16
<p>
Expanding hypergeometric functions in polynomial closed form before evaluating them can lead to numerical cancellation. Compare for example
</p>
<pre class="wiki">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
</pre><p>
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.)
</p>
TicketrwsTue, 20 Jan 2015 07:21:41 GMT
https://trac.sagemath.org/ticket/17066#comment:17
https://trac.sagemath.org/ticket/17066#comment:17
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/17066#comment:16" title="Comment 16">mmezzarobba</a>:
</p>
<blockquote class="citation">
<p>
Won't that be a problem with this ticket?
</p>
</blockquote>
<p>
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
</p>
<p>
EDIT: nonsense about <code>n()</code> removed
</p>
Ticketfredrik.johanssonTue, 20 Jan 2015 12:22:12 GMT
https://trac.sagemath.org/ticket/17066#comment:18
https://trac.sagemath.org/ticket/17066#comment:18
<p>
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.
</p>
TicketmmezzarobbaTue, 20 Jan 2015 12:54:21 GMT
https://trac.sagemath.org/ticket/17066#comment:19
https://trac.sagemath.org/ticket/17066#comment:19
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/17066#comment:17" title="Comment 17">rws</a>:
</p>
<blockquote class="citation">
<p>
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
</p>
</blockquote>
<p>
If (as you seem to say in your post to <code>sage-devel</code>) the result of the ticket is that
</p>
<pre class="wiki">sage: hypergeometric([1,-100],[3],x)
</pre><p>
returns a polynomial while
</p>
<pre class="wiki">sage: hypergeometric([1,-100],[3],0.7)
</pre><p>
still calls a numerically stable evaluation algorithm, then I don't think there is a problem. (Ideally, <code> hypergeometric([1,-100],[3],x,hold=True).subs(x=0.7)</code> should also call the numerically stable evaluation algorithm, though.)
</p>
<p>
What I was concerned about was that the implementation of <code>n()</code> for an "atomic-looking" expression with dedicated numerical evaluation code would end up using a numerically unstable formula <em>internally</em> 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.
</p>
TicketrwsTue, 20 Jan 2015 13:38:27 GMT
https://trac.sagemath.org/ticket/17066#comment:20
https://trac.sagemath.org/ticket/17066#comment:20
<p>
It's not that I don't see your point. The real fix for this cannot be provided by <code>hypergeometric</code> though but must be done in <code>subs</code>, if at all.
</p>
Ticketfredrik.johanssonTue, 20 Jan 2015 13:46:04 GMT
https://trac.sagemath.org/ticket/17066#comment:21
https://trac.sagemath.org/ticket/17066#comment:21
<p>
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.
</p>
TicketgitTue, 20 Jan 2015 17:13:00 GMTcommit changed
https://trac.sagemath.org/ticket/17066#comment:22
https://trac.sagemath.org/ticket/17066#comment:22
<ul>
<li><strong>commit</strong>
changed from <em>55daf81080f46c54429e629b65e8049acfdb6d00</em> to <em>0e0fef65a021a9b0a9db30a2bc6d300ede68bf14</em>
</li>
</ul>
<p>
Branch pushed to git repo; I updated commit sha1. New commits:
</p>
<table class="wiki">
<tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=290b47f1c5f1ca7c8afcf6985e1942f563338671"><span class="icon"></span>290b47f</a></td><td><code>Merge branch 'develop' into t/17066/work_around_mpmath_problem_with_hypergeometric___zeroes</code>
</td></tr><tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=0e0fef65a021a9b0a9db30a2bc6d300ede68bf14"><span class="icon"></span>0e0fef6</a></td><td><code>17066: add warning</code>
</td></tr></table>
TicketkcrismanTue, 20 Jan 2015 17:42:29 GMT
https://trac.sagemath.org/ticket/17066#comment:23
https://trac.sagemath.org/ticket/17066#comment:23
<blockquote class="citation">
<p>
Expanding hypergeometric functions in polynomial closed form before evaluating them can lead to numerical cancellation. Compare for example
</p>
<pre class="wiki">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
</pre><p>
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.)
</p>
</blockquote>
<p>
Just so I follow, what you are saying is that the polynomial (rational function, I guess) given will be <em>correct</em>, 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 <em>any</em> 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.
</p>
<blockquote class="citation">
<p>
The real fix for this cannot be provided by hypergeometric though but must be done in subs, if at all.
</p>
</blockquote>
<p>
Well, it could provide a good reason <em>not</em> 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?
</p>
TicketmmezzarobbaTue, 20 Jan 2015 18:06:56 GMT
https://trac.sagemath.org/ticket/17066#comment:24
https://trac.sagemath.org/ticket/17066#comment:24
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/17066#comment:23" title="Comment 23">kcrisman</a>:
</p>
<blockquote class="citation">
<p>
Just so I follow, what you are saying is that the polynomial (rational function, I guess)
</p>
</blockquote>
<p>
No, it is really a polynomial—but one of degree 100, with large coefficients that alternate in sign.
</p>
<blockquote class="citation">
<p>
given will be <em>correct</em>, but that substituting floating point numbers are likely to be inaccurate unless carried out to a very high degree of precision?
</p>
</blockquote>
<p>
Precisely. (Well, not necessarily a <em>very</em> high precision, but probably something like 80-100 guard bits in my example.)
</p>
TicketkcrismanTue, 20 Jan 2015 18:12:55 GMT
https://trac.sagemath.org/ticket/17066#comment:25
https://trac.sagemath.org/ticket/17066#comment:25
<blockquote class="citation">
<blockquote class="citation">
<p>
Just so I follow, what you are saying is that the polynomial (rational function, I guess)
</p>
</blockquote>
<p>
No, it is really a polynomial—but one of degree 100, with large coefficients that alternate in sign.
</p>
</blockquote>
<p>
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.
</p>
<blockquote class="citation">
<blockquote class="citation">
<p>
given will be <em>correct</em>, but that substituting floating point numbers are likely to be inaccurate unless carried out to a very high degree of precision?
</p>
</blockquote>
<p>
Precisely. (Well, not necessarily a <em>very</em> high precision, but probably something like 80-100 guard bits in my example.)
</p>
</blockquote>
<p>
Anyway, 53 default yields madness, as we see in <a class="ticket" href="https://trac.sagemath.org/ticket/17066#comment:16" title="Comment 16">comment:16</a>. Not good, but I'm not sure what the "right" answer is; I'm leaning toward saying we should <em>not</em> 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?
</p>
TicketrwsTue, 24 Feb 2015 09:45:57 GMTstatus changed
https://trac.sagemath.org/ticket/17066#comment:26
https://trac.sagemath.org/ticket/17066#comment:26
<ul>
<li><strong>status</strong>
changed from <em>needs_review</em> to <em>needs_work</em>
</li>
</ul>
<pre class="wiki">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
</pre>
Ticket