Sage: Ticket #10268: adding GiNaC method to simplify_rational
https://trac.sagemath.org/ticket/10268
<p>
Currently simplify_rational() only offers 3 Maxima methods. GiNaC offers another possibility via its normal() method. This issue is discussed here
</p>
<p>
<a class="ext-link" href="http://groups.google.com/group/sage-devel/browse_thread/thread/843c17dcbd9c2958"><span class="icon"></span>http://groups.google.com/group/sage-devel/browse_thread/thread/843c17dcbd9c2958</a>
</p>
<p>
I have a patch and a benchmark but need to redownload sage because I am getting unrelated doctest failures with or without the patch.
</p>
<p>
EDIT: All tests pass now with the attached patch, as they should because the default behavior is not changed. Also, I am attaching a benchmark script using random rational expressions that simplify to 1. In this benchmark, the GiNaC option is about 10 times faster than the default option (Maxima's fullratsimp, without utilizing libraryness).
</p>
<p>
One limitation of this patch is that it does not support Maxima's map option. GiNaC has a map function, but utilizing it from sage would require a bit more effort.
</p>
en-usSagehttps://trac.sagemath.org/chrome/site/logo_sagemath_trac.png
https://trac.sagemath.org/ticket/10268
Trac 1.1.6bgoodriMon, 15 Nov 2010 15:48:09 GMTattachment set
https://trac.sagemath.org/ticket/10268
https://trac.sagemath.org/ticket/10268
<ul>
<li><strong>attachment</strong>
set to <em>trac_10268_enhance_simplify_rational.patch</em>
</li>
</ul>
TicketbgoodriMon, 15 Nov 2010 15:52:28 GMTdescription changed
https://trac.sagemath.org/ticket/10268#comment:1
https://trac.sagemath.org/ticket/10268#comment:1
<ul>
<li><strong>description</strong>
modified (<a href="/ticket/10268?action=diff&version=1">diff</a>)
</li>
</ul>
TicketbgoodriMon, 15 Nov 2010 22:24:42 GMTattachment set
https://trac.sagemath.org/ticket/10268
https://trac.sagemath.org/ticket/10268
<ul>
<li><strong>attachment</strong>
set to <em>test.ginsh</em>
</li>
</ul>
TicketbgoodriMon, 15 Nov 2010 22:24:55 GMTattachment set
https://trac.sagemath.org/ticket/10268
https://trac.sagemath.org/ticket/10268
<ul>
<li><strong>attachment</strong>
set to <em>test.sage</em>
</li>
</ul>
TicketbgoodriMon, 15 Nov 2010 22:27:06 GMT
https://trac.sagemath.org/ticket/10268#comment:2
https://trac.sagemath.org/ticket/10268#comment:2
<p>
Okay, actually my patch does not work so well for the example that motivated it. I am now attaching a test.ginsh file that defines a rational expression, substitutes a variable with another big rational expression, calls normal and quits. The GiNaC shell finishes in about 1 minute on my laptop (in a shell execute: "time ginsh test.ginsh"). But after applying my patch to sage and then trying to do the equivalent thing via test.sage, it can go for hours without finishing. So, something is badly wrong, possibly my patch.
</p>
TicketkcrismanTue, 16 Nov 2010 01:57:43 GMT
https://trac.sagemath.org/ticket/10268#comment:3
https://trac.sagemath.org/ticket/10268#comment:3
<p>
I'm not a Cython expert, but maybe should you only use the <code>sig_on</code> business inside the else? It would make sense to start off with method='normal' and then do Ginac, otherwise do Maxima; there is no need to create the <code>GEx</code> unless you are actually going to use it, I would think. I don't know if this would have anything to do with the hang, but it's worth a shot.
</p>
TicketbgoodriTue, 16 Nov 2010 04:33:06 GMT
https://trac.sagemath.org/ticket/10268#comment:4
https://trac.sagemath.org/ticket/10268#comment:4
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/10268#comment:3" title="Comment 3">kcrisman</a>:
</p>
<p>
I am just learning sage, but it seems that the compiler does not like the GEx to be declared inside a conditional statement, which makes sense. The _sig_on and _sig_off thing I think is for catching segfaults, which doesn't seem to be a problem and when I comment those out, the behavior is the same. A slight possibility is the fact that when I use the GiNaC shell directly it is the most recent version, whereas Pynac forked off an older version, but the normal function has been in GiNaC for a long, long time.
</p>
<p>
More interesting is that when I interrupt sage, I get this traceback
</p>
<pre class="wiki">
KeyboardInterrupt Traceback (most recent call last)
/media/disk30/sage-4.6/<ipython console> in <module>()
/media/disk30/sage-4.6/local/lib/python2.6/site-packages/sage/symbolic/expression.so in sage.symbolic.expression.Expression.simplify_rational (sage/symbolic/expression.cpp:23989)()
/media/disk30/sage-4.6/local/lib/python2.6/site-packages/sage/symbolic/pynac.so in sage.symbolic.pynac.py_gcd (sage/symbolic/pynac.cpp:6440)()
/media/disk30/sage-4.6/local/lib/python2.6/site-packages/sage/rings/arith.pyc in gcd(a, b, **kwargs)
1363 sigma = Sigma()
1364
-> 1365 def gcd(a, b=None, **kwargs):
1366 r"""
1367 The greatest common divisor of a and b, or if a is a list and b is
/media/disk30/sage-4.6/local/lib/python2.6/site-packages/sage/interfaces/get_sigs.pyc in my_sigint(x, n)
7
8 def my_sigint(x, n):
----> 9 raise KeyboardInterrupt
10
11 def my_sigfpe(x, n):
KeyboardInterrupt:
</pre><p>
Why would it fall into the gcd function from /media/disk30/sage-4.6/local/lib/python2.6/site-packages/sage/rings/arith.pyc? The patch does not call it directly, and it is a waste because normal in GiNaC already cancels the greatest common factor from the numerator and the denominator. And then a related question is why does gcd seem to hang?
</p>
TicketbgoodriTue, 16 Nov 2010 16:12:06 GMTattachment set
https://trac.sagemath.org/ticket/10268
https://trac.sagemath.org/ticket/10268
<ul>
<li><strong>attachment</strong>
set to <em>bench.sage</em>
</li>
</ul>
TicketburcinTue, 16 Nov 2010 16:12:45 GMTstatus changed; author set
https://trac.sagemath.org/ticket/10268#comment:5
https://trac.sagemath.org/ticket/10268#comment:5
<ul>
<li><strong>status</strong>
changed from <em>new</em> to <em>needs_work</em>
</li>
<li><strong>author</strong>
set to <em>Ben Goodrich</em>
</li>
</ul>
<p>
I'm really happy to see some effort to use pynac/ginac to replace functionality we normally use maxima for. Unfortunately this is a really busy period for me so I can't help much. Thanks a lot for your effort Ben.
</p>
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/10268#comment:4" title="Comment 4">bgoodri</a>:
</p>
<blockquote class="citation">
<p>
Why would it fall into the gcd function from /media/disk30/sage-4.6/local/lib/python2.6/site-packages/sage/rings/arith.pyc? The patch does not call it directly, and it is a waste because normal in GiNaC already cancels the greatest common factor from the numerator and the denominator. And then a related question is why does gcd seem to hang?
</p>
</blockquote>
<p>
The <code>numeric::gcd()</code> method calls <code>sage.symbolic.pynac.py_gcd()</code>. See here:
</p>
<p>
<a class="ext-link" href="http://pynac.sagemath.org/hg/file/b233d9dadcfa/ginac/numeric.cpp#l2526"><span class="icon"></span>http://pynac.sagemath.org/hg/file/b233d9dadcfa/ginac/numeric.cpp#l2526</a>
</p>
<p>
It could be that our gcd() function doesn't work exactly like CLN's <code>gcd()</code> function which is used originally by ginac. This would effect the termination criteria used in the multivariate gcd code in ginac/pynac.
</p>
<p>
I haven't looked into the functionality in <code>normal.cpp</code> much, but one of William's goals was to make it call Singular (or the Factory library) to factor multivariate polynomials instead of the code in ginac. This library generally performs much better and it is actively being developed.
</p>
<hr />
<p>
BTW, kcrisman was right about his comment on the use of <code>sig_on/sig_off</code>. These functions are used to override the signal handlers so we can catch <code>CTRL-C</code> in long running library code. You don't need them around the declaration of <code>GEX</code>, but you should use them in the call to <code>normal()</code>. AFAIR, Jeroen Demeyer recently wrote a section on this in the developers guide, but I don't have a link handy right now. The new call signature for these might be <code>sig_on()</code>, like a function call.
</p>
TicketbgoodriTue, 16 Nov 2010 17:28:38 GMTcc set
https://trac.sagemath.org/ticket/10268#comment:6
https://trac.sagemath.org/ticket/10268#comment:6
<ul>
<li><strong>cc</strong>
<em>was</em> added
</li>
</ul>
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/10268#comment:5" title="Comment 5">burcin</a>: cc'ing William for clarification
</p>
<blockquote class="citation">
<p>
I haven't looked into the functionality in <code>normal.cpp</code> much, but one of William's goals was to make it call Singular (or the Factory library) to factor multivariate polynomials instead of the code in ginac. This library generally performs much better and it is actively being developed.
</p>
</blockquote>
<p>
Do you really want to do full factorization in <code>simplify_rational()</code>. I think neither Maxima nor GiNaC do that, only square-free factorization and gcd cancellation. We could add an option to do full factorization of the numerator and denominator before returning. If so, would it make more sense to first backport the functionality in factor.cpp from GiNaC 1.5.x to the pynac fork than to code a pynac-libSingular link?
</p>
<p>
Also, I think this might be a bit separate from the issue I was hitting. When I ran test.sage last night under trace() with the enter key wedged down, by the morning it had called gcd() over 30,000 times and hadn't even passed the rational expression to GiNaC yet. This is a waste because GiNaC's normal() function was going to do 1 gcd cancellation anyway. So, it seems what we need is an option to prevent sage from trying to find the gcd of every subexpression.
</p>
TicketbgoodriTue, 16 Nov 2010 17:38:14 GMT
https://trac.sagemath.org/ticket/10268#comment:7
https://trac.sagemath.org/ticket/10268#comment:7
<p>
The original bench.sage was not very appropriate because sage was simplifying the rational expression to 1 before passing it to Maxima or GiNaC. So the difference in speed primarily reflected the difference between interacting via pexpect and interacting via a library. The revised bench.sage avoids this and there is about a 7x speedup. However, sage is repeatedly calling gcd() automatically, and the performance would probably jump if we could avoid that somehow.
</p>
TicketbgoodriWed, 17 Nov 2010 02:25:01 GMT
https://trac.sagemath.org/ticket/10268#comment:8
https://trac.sagemath.org/ticket/10268#comment:8
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/10268#comment:5" title="Comment 5">burcin</a>:
</p>
<blockquote class="citation">
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/10268#comment:4" title="Comment 4">bgoodri</a>:
</p>
<blockquote class="citation">
<p>
And then a related question is why does gcd seem to hang?
</p>
</blockquote>
<p>
The <code>numeric::gcd()</code> method calls <code>sage.symbolic.pynac.py_gcd()</code>. See here: <a class="ext-link" href="http://pynac.sagemath.org/hg/file/b233d9dadcfa/ginac/numeric.cpp#l2526"><span class="icon"></span>http://pynac.sagemath.org/hg/file/b233d9dadcfa/ginac/numeric.cpp#l2526</a> It could be that our gcd() function doesn't work exactly like CLN's <code>gcd()</code> function which is used originally by ginac. This would effect the termination criteria used in the multivariate gcd code in ginac/pynac. <br />
</p>
</blockquote>
<p>
Okay, I've been contributing to the confusion. Now I see what you meant: When sage calls Pynac's normal() function, Pynac calls "its" gcd() function, which is actually sage's gcd() function. So, all the calls to gcd() are expected behavior, and the question becomes why doesn't Pynac get them over with and terminate in 30 seconds like GiNaC does with the CLN implementation of gcd()?
</p>
<p>
For reference when I switch on the statistics bookkeeping in the (latest) GiNaC source:
</p>
<pre class="wiki">goodrich@Y560:/tmp$ time ginsh test.ginsh
ginsh - GiNaC Interactive Shell (ginac V1.5.8)
__, _______ Copyright (C) 1999-2010 Johannes Gutenberg University Mainz,
(__) * | Germany. This is free software with ABSOLUTELY NO WARRANTY.
._) i N a C | You are welcome to redistribute it under certain conditions.
<-------------' For details type `warranty;'.
Type ?? for a list of help topics.
gcd() called 56331 times
sr_gcd() called 0 times
heur_gcd() called 1612 times
heur_gcd() failed 1 times
real 0m34.025s
user 0m33.245s
sys 0m0.602s
</pre><p>
Okay, that is what is supposed to happen. With the current Pynac, I interrupt after an hour and
</p>
<pre class="wiki">sage: quit;
Exiting Sage (CPU time 62m29.29s, Wall time 64m3.75s).
gcd() called 57537 times
sr_gcd() called 19 times
heur_gcd() called 1936 times
heur_gcd() failed 19 times
</pre><p>
So, it looks as if Pynac is hanging at or toward the end, and it experiences many more failures in the heur_gcd() routine. I guess I should be looking at the gcd heuristics then. Any ideas come to mind?
</p>
TicketbgoodriWed, 17 Nov 2010 08:26:27 GMT
https://trac.sagemath.org/ticket/10268#comment:9
https://trac.sagemath.org/ticket/10268#comment:9
<p>
Bah, it was a bug in GiNaC that was fixed by this recent GiNaC commit
</p>
<p>
<a class="ext-link" href="http://www.ginac.de/ginac.git?p=ginac.git;a=commit;h=edf1ae46a926d0a718063c149b78c1b9a7ec2043"><span class="icon"></span>http://www.ginac.de/ginac.git?p=ginac.git;a=commit;h=edf1ae46a926d0a718063c149b78c1b9a7ec2043</a>
</p>
<p>
I can bring the bug back by reverting it.
</p>
<p>
However, the patch touches code that was only added to the 1.5.x branch of GiNaC, so we can't just apply that patch to the pynac spkg. I guess the logic that this patch fixes was also wrong somewhere in the 1.4.x branch of GiNaC. But I'm too tired and frustrated to look into it right now.
</p>
TicketburcinThu, 18 Nov 2010 15:47:49 GMT
https://trac.sagemath.org/ticket/10268#comment:10
https://trac.sagemath.org/ticket/10268#comment:10
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/10268#comment:6" title="Comment 6">bgoodri</a>:
</p>
<blockquote class="citation">
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/10268#comment:5" title="Comment 5">burcin</a>: cc'ing William for clarification
</p>
<blockquote class="citation">
<p>
I haven't looked into the functionality in <code>normal.cpp</code> much, but one of William's goals was to make it call Singular (or the Factory library) to factor multivariate polynomials instead of the code in ginac. This library generally performs much better and it is actively being developed.
</p>
</blockquote>
<p>
Do you really want to do full factorization in <code>simplify_rational()</code>. I think neither Maxima nor GiNaC do that, only square-free factorization and gcd cancellation. We could add an option to do full factorization of the numerator and denominator before returning. If so, would it make more sense to first backport the functionality in factor.cpp from GiNaC 1.5.x to the pynac fork than to code a pynac-libSingular link?
</p>
</blockquote>
<p>
Sorry for the confusion. I meant to say gcd.
</p>
<p>
Now that there is a separate ticket for the bug in the pynac gcd, <a class="closed ticket" href="https://trac.sagemath.org/ticket/10284" title="defect: Infinite loop in gcd() via pynac-0.2.1 (closed: invalid)">#10284</a>, I will post my response to the other questions there.
</p>
TicketjdemeyerTue, 13 Aug 2013 15:35:53 GMTmilestone changed
https://trac.sagemath.org/ticket/10268#comment:11
https://trac.sagemath.org/ticket/10268#comment:11
<ul>
<li><strong>milestone</strong>
changed from <em>sage-5.11</em> to <em>sage-5.12</em>
</li>
</ul>
TicketeviatarbachMon, 26 Aug 2013 17:28:47 GMTcc changed
https://trac.sagemath.org/ticket/10268#comment:12
https://trac.sagemath.org/ticket/10268#comment:12
<ul>
<li><strong>cc</strong>
<em>eviatarbach</em> added
</li>
</ul>
Ticketvbraun_spamThu, 30 Jan 2014 21:20:52 GMTmilestone changed
https://trac.sagemath.org/ticket/10268#comment:13
https://trac.sagemath.org/ticket/10268#comment:13
<ul>
<li><strong>milestone</strong>
changed from <em>sage-6.1</em> to <em>sage-6.2</em>
</li>
</ul>
Ticketvbraun_spamTue, 06 May 2014 15:20:58 GMTmilestone changed
https://trac.sagemath.org/ticket/10268#comment:14
https://trac.sagemath.org/ticket/10268#comment:14
<ul>
<li><strong>milestone</strong>
changed from <em>sage-6.2</em> to <em>sage-6.3</em>
</li>
</ul>
Ticketvbraun_spamSun, 10 Aug 2014 16:51:03 GMTmilestone changed
https://trac.sagemath.org/ticket/10268#comment:15
https://trac.sagemath.org/ticket/10268#comment:15
<ul>
<li><strong>milestone</strong>
changed from <em>sage-6.3</em> to <em>sage-6.4</em>
</li>
</ul>
TicketrwsSat, 11 Jun 2016 13:42:26 GMT
https://trac.sagemath.org/ticket/10268#comment:16
https://trac.sagemath.org/ticket/10268#comment:16
<p>
Pynac's <code>normal()</code> was made accessible by <code>Expression</code> in <a class="closed ticket" href="https://trac.sagemath.org/ticket/12068" title="enhancement: Numerator for symbolic expression shouldn't use maxima (closed: fixed)">#12068</a>, so this patch can take advantage of it. The gcd error can be worked around optionally with Pynac-0.6.7 and libgiac if it is installed (<a class="closed ticket" href="https://trac.sagemath.org/ticket/20742" title="defect: Upgrade to pynac-0.6.7 (closed: fixed)">#20742</a>). This speeds up the GCD by an order of mag versus Pynac and two orders of mag versus pexpect-Maxima.
</p>
Ticket