Sage: Ticket #16964: Speed up comparisons in QQbar
https://trac.sagemath.org/ticket/16964
<p>
When computing the variety of some ideal, then excessive amounts of time are apparently spent sorting the solutions. (In one example, the variety was essentially computed in 7½ hours but the sorting wasn't finished even 1½ hours after that.) This is due to the fact that comparison between complex algebraic numbers is done lexicographically, which means that quite often the real parts of complex conjugate numbers will have to be compared for equality, which can be a very costly operation.
</p>
<p>
Originally the computation even resulted in a Pari exception (“not enough precomputed primes”). This no longer occurs (since the pari upgrade from <a class="closed ticket" href="https://trac.sagemath.org/ticket/15767" title="#15767: enhancement: Upgrade PARI to 2.7.1 (closed: fixed)">#15767</a>). So the focus of this ticket is now the excessive amount of time required for comparisons, even without an exception.
</p>
en-usSagehttps://trac.sagemath.org/chrome/site/logo_sagemath_trac.png
https://trac.sagemath.org/ticket/16964
Trac 1.2Jeroen DemeyerThu, 11 Sep 2014 07:00:33 GMT
https://trac.sagemath.org/ticket/16964#comment:1
https://trac.sagemath.org/ticket/16964#comment:1
<p>
Please mention the version of Sage, in particular whether or not <a class="closed ticket" href="https://trac.sagemath.org/ticket/15767" title="#15767: enhancement: Upgrade PARI to 2.7.1 (closed: fixed)">#15767</a> was applied.
</p>
TicketMartin von GagernFri, 12 Sep 2014 07:00:24 GMT
https://trac.sagemath.org/ticket/16964#comment:2
https://trac.sagemath.org/ticket/16964#comment:2
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/16964#comment:1" title="Comment 1">jdemeyer</a>:
</p>
<blockquote class="citation">
<p>
Please mention the version of Sage, in particular whether or not <a class="closed ticket" href="https://trac.sagemath.org/ticket/15767" title="#15767: enhancement: Upgrade PARI to 2.7.1 (closed: fixed)">#15767</a> was applied.
</p>
</blockquote>
<p>
This was sage 6.3 on Gentoo. With current develop (6.4.beta2) at least it doesn't fail as “quickly”: where the previous computation was 8 hours all in all up to the exception, my patched version without sorting took 7½ hours to compute the list, and after that sage has now been busy sorting that list for 1½ hours and still isn't done. I'll report back when it is, or when it gave up, but I'd say some need for action remains, even if the exception no longer occurs in this form.
</p>
TicketMartin von GagernFri, 12 Sep 2014 07:22:58 GMT
https://trac.sagemath.org/ticket/16964#comment:3
https://trac.sagemath.org/ticket/16964#comment:3
<p>
Here is a reproducing example which at least demonstrates that comparisons take <em>way</em> longer than they should:
</p>
<pre class="wiki">sage: r = QQ[x](69721504*x^8 + 251777664*x^6 + 329532012*x^4 + 184429548*x^2 + 37344321).roots(QQbar, False)
sage: r
[-0.0221204634374360? - 1.090991904211621?*I,
-0.0221204634374360? + 1.090991904211621?*I,
-0.8088604911480535?*I,
-0.7598602580415435?*I,
0.7598602580415435?*I,
0.8088604911480535?*I,
0.0221204634374360? - 1.090991904211621?*I,
0.0221204634374360? + 1.090991904211621?*I]
sage: [r[0], r[1]].sort()
</pre><p>
This is because the comparison of the real parts takes like forever. Which in turn is because the computation of its minimal polynomial takes forever.
</p>
<p>
Looking at the set of all zeros, I can see that there are 4 clearly distinct real parts, and each comes with a pair of conjugate solutions since the polynomial has real coefficients. This is enough to conclude that if the intervals for two real parts overlap, then they must be equal and I don't have to do an exact computation for this. Should we try to implement some of this reasoning as a special case for <code>AlgebraicNumber.__cmp__</code>, for the case where the descriptor is exact and the minpoly is the same?
</p>
TicketJeroen DemeyerSun, 14 Sep 2014 17:43:08 GMTsummary changed
https://trac.sagemath.org/ticket/16964#comment:4
https://trac.sagemath.org/ticket/16964#comment:4
<ul>
<li><strong>summary</strong>
changed from <em>MPolynomialIdeal_singular_repr.variety: sage.libs.pari.gen.PariError: not enough precomputed primes</em> to <em>Fix comparisons in QQbar</em>
</li>
</ul>
TicketMartin von GagernMon, 06 Oct 2014 08:27:36 GMT
https://trac.sagemath.org/ticket/16964#comment:5
https://trac.sagemath.org/ticket/16964#comment:5
<p>
Looking for ways to address this, I noticed that the a lot of time apparently is spent inside
</p>
<pre class="wiki">class ANBinaryExpr(ANDescr):
⋮
def exactify(self):
⋮
gen = left._exact_field().union(right._exact_field())
⋮
</pre><p>
I wonder whether we can avoid that union for the case where both fields have the same defining polynomial. I wonder whether we can assume that the root of the field will form a <a class="ext-link" href="http://en.wikipedia.org/wiki/Algebraic_number_field#Power_basis"><span class="icon"></span>power basis</a>, and if so, whether there is any reasonably cheap way to find the conversion between different power bases, so we could express one root in terms of another.
</p>
<p>
Since I don't have any good ideas how to achieve this, my best idea still is tacking this at the <code>__cmp__</code> level, but if anyone has an idea for solving this more generic issue, that would be really great since it would help other computations as well. So I'm sharing my thoughts.
</p>
<p>
Here is what I've tried and discarded, so you can avoid that same cul de sac. I started by writing down a generic linear combination w = a₀ + a₁z + a₂z² + … and then computed p(w) reduced by p(z), where p is the polynomial of the field. This gives a polynomial in z, and if enough powers of z are irrational then all the coefficients have to be zero if w is a root and the a_i are to be rational. So this gave me d conditions on these a₀ through a_{d-1}, which I could combine into an ideal and try to compute a variety. But that variety computation takes like forever in the above example, so this generic approach of finding other roots is not feasible in this fashion. Been there, tried that and discarded it.
</p>
TicketJeroen DemeyerSat, 13 Dec 2014 22:19:12 GMTpriority, component, summary changed
https://trac.sagemath.org/ticket/16964#comment:6
https://trac.sagemath.org/ticket/16964#comment:6
<ul>
<li><strong>priority</strong>
changed from <em>major</em> to <em>critical</em>
</li>
<li><strong>component</strong>
changed from <em>algebraic geometry</em> to <em>number fields</em>
</li>
<li><strong>summary</strong>
changed from <em>Fix comparisons in QQbar</em> to <em>Speed up comparisons in QQbar</em>
</li>
</ul>
TicketJeroen DemeyerSat, 13 Dec 2014 22:20:12 GMTdescription, milestone changed
https://trac.sagemath.org/ticket/16964#comment:7
https://trac.sagemath.org/ticket/16964#comment:7
<ul>
<li><strong>description</strong>
modified (<a href="/ticket/16964?action=diff&version=7">diff</a>)
</li>
<li><strong>milestone</strong>
changed from <em>sage-6.4</em> to <em>sage-6.5</em>
</li>
</ul>
TicketMartin von GagernMon, 15 Dec 2014 10:39:04 GMTchangetime, time changed; branch set
https://trac.sagemath.org/ticket/16964#comment:8
https://trac.sagemath.org/ticket/16964#comment:8
<ul>
<li><strong>changetime</strong>
changed from <em>Dec 13, 2014, 10:20:12 PM</em> to <em>Dec 13, 2014, 10:20:12 PM</em>
</li>
<li><strong>branch</strong>
set to <em>u/gagern/ticket/16964</em>
</li>
<li><strong>time</strong>
changed from <em>Sep 11, 2014, 6:23:53 AM</em> to <em>Sep 11, 2014, 6:23:53 AM</em>
</li>
</ul>
TicketMartin von GagernMon, 15 Dec 2014 10:47:01 GMTstatus changed; commit, author set
https://trac.sagemath.org/ticket/16964#comment:9
https://trac.sagemath.org/ticket/16964#comment:9
<ul>
<li><strong>status</strong>
changed from <em>new</em> to <em>needs_review</em>
</li>
<li><strong>commit</strong>
set to <em>9f6fdc209d62e055df4bdd38f9f946332fbf3842</em>
</li>
<li><strong>author</strong>
set to <em>Martin von Gagern</em>
</li>
</ul>
<p>
Here is my first stab at the approach outlined in <a class="ticket" href="https://trac.sagemath.org/ticket/16964#comment:3" title="Comment 3">comment:3</a>. If you have a nicer way to handle things, feel free to suggest an alternative. If you think that my special case should not call <code>minpoly()</code> but instead examine whether the descriptor is <code>ANRoot</code> with matching polynomial, please state your reason for this argument. Likewise, if you think I should also handle more than one real value or conjugate pair for a given real component, then I'd love to hear how you'd implement this without making the code too hard to read and maintain.
</p>
<p>
As it is, I consider my code not particularly beautiful, and perhaps not the best solution either, but it is <em>way</em> better than the current state of affairs, and since we have other things depending on this, e.g. <a class="closed ticket" href="https://trac.sagemath.org/ticket/14239" title="#14239: enhancement: symbolic radical expression for algebraic number (closed: fixed)">#14239</a>, I'd like to see this merged pretty soon and possible improvements dealt with in a later ticket.
</p>
<hr />
<p>
New commits:
</p>
<table class="wiki">
<tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=9f6fdc209d62e055df4bdd38f9f946332fbf3842"><span class="icon"></span>9f6fdc2</a></td><td><code>Faster qqbar comparison for case with equal minpoly.</code>
</td></tr></table>
TicketgitMon, 15 Dec 2014 11:04:49 GMTcommit changed
https://trac.sagemath.org/ticket/16964#comment:10
https://trac.sagemath.org/ticket/16964#comment:10
<ul>
<li><strong>commit</strong>
changed from <em>9f6fdc209d62e055df4bdd38f9f946332fbf3842</em> to <em>2c4ac1be555a9c81dee84acfc7d7c2cb2602db22</em>
</li>
</ul>
<p>
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
</p>
<table class="wiki">
<tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=2c4ac1be555a9c81dee84acfc7d7c2cb2602db22"><span class="icon"></span>2c4ac1b</a></td><td><code>Faster qqbar comparison for case with equal minpoly.</code>
</td></tr></table>
TicketJeroen DemeyerMon, 15 Dec 2014 22:11:44 GMT
https://trac.sagemath.org/ticket/16964#comment:11
https://trac.sagemath.org/ticket/16964#comment:11
<p>
Just a quick comment: I like your general approach, but I have to check the details...
</p>
TicketgitThu, 08 Jan 2015 19:47:19 GMTcommit changed
https://trac.sagemath.org/ticket/16964#comment:12
https://trac.sagemath.org/ticket/16964#comment:12
<ul>
<li><strong>commit</strong>
changed from <em>2c4ac1be555a9c81dee84acfc7d7c2cb2602db22</em> to <em>e7acd96040fb43eebed0e3df0f1b1c7613d49bc1</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=e7acd96040fb43eebed0e3df0f1b1c7613d49bc1"><span class="icon"></span>e7acd96</a></td><td><code>Fix term in documentation.</code>
</td></tr></table>
TicketMartin von GagernMon, 26 Jan 2015 09:36:43 GMTdescription changed
https://trac.sagemath.org/ticket/16964#comment:13
https://trac.sagemath.org/ticket/16964#comment:13
<ul>
<li><strong>description</strong>
modified (<a href="/ticket/16964?action=diff&version=13">diff</a>)
</li>
</ul>
<p>
Changing description to turn focus from pari exception to performance.
</p>
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/16964#comment:11" title="Comment 11">jdemeyer</a>:
</p>
<blockquote class="citation">
<p>
I like your general approach, but I have to check the details...
</p>
</blockquote>
<p>
Have you found time for a closer look? Since this is currently the only ticket in the review queue with priority critical, perhaps we should try to get this into Sage 6.5?
</p>
TicketVincent DelecroixSat, 28 Feb 2015 14:15:44 GMTstatus changed
https://trac.sagemath.org/ticket/16964#comment:14
https://trac.sagemath.org/ticket/16964#comment:14
<ul>
<li><strong>status</strong>
changed from <em>needs_review</em> to <em>needs_info</em>
</li>
</ul>
<p>
Hello,
</p>
<p>
This improvement is really cool!
</p>
<p>
EDIT: the thing below is basically one part of <a class="ticket" href="https://trac.sagemath.org/ticket/16964#comment:9" title="Comment 9">comment:9</a>...
</p>
<p>
Couldn't we do a little bit better? In some cases we might identify numerically pairs of conjugated roots (and possibly one real root) and still be able to decide which one is which. There exists polynomials whose roots have the same real part:
</p>
<pre class="wiki">sage: x = polygen(ZZ)
sage: P = x^4 - 4*x^3 + 9*x^2 - 10*x + 5
sage: P.roots(CC,False)
[1.00000000000000 - 1.61803398874989*I,
1.00000000000000 - 0.618033988749895*I,
1.00000000000000 + 0.618033988749895*I,
1.00000000000000 + 1.61803398874989*I]
</pre><p>
With the implementation proposed in this ticket, the comparison of the roots of the above polynomial falls back to the generic implementation. We could at least compare the conjugated ones without exactification of the real part. But we can leave that for a further ticket.
</p>
<p>
Vincent
</p>
TicketVincent DelecroixSat, 28 Feb 2015 14:34:00 GMT
https://trac.sagemath.org/ticket/16964#comment:15
https://trac.sagemath.org/ticket/16964#comment:15
<p>
And you should also avoid this special case if <code>-self._value.imag()</code> and <code>other._value.imag()</code> do not overlap. That would potentially save a call to <code>minpoly</code>.
</p>
TicketVincent DelecroixSat, 28 Feb 2015 14:48:55 GMTstatus, commit, branch changed
https://trac.sagemath.org/ticket/16964#comment:16
https://trac.sagemath.org/ticket/16964#comment:16
<ul>
<li><strong>status</strong>
changed from <em>needs_info</em> to <em>needs_review</em>
</li>
<li><strong>commit</strong>
changed from <em>e7acd96040fb43eebed0e3df0f1b1c7613d49bc1</em> to <em>1294bf20bb6fe915a5371318cb8f9feacfa3f86a</em>
</li>
<li><strong>branch</strong>
changed from <em>u/gagern/ticket/16964</em> to <em>u/vdelecroix/16964</em>
</li>
</ul>
<p>
I pushed a review commit. Tell me what you think. It still lacks a nice complicated example... (I am looking for it)
</p>
<hr />
<p>
New commits:
</p>
<table class="wiki">
<tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=8c0cf63c0e86d0cfc042149901ef9cfcde4d14ac"><span class="icon"></span>8c0cf63</a></td><td><code>trac #17863: remove stuff from src/ext</code>
</td></tr><tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=2777e6bb12d802a64528bc37f85fa7be561c5b29"><span class="icon"></span>2777e6b</a></td><td><code>Merge branch 'u/gagern/ticket/16964' of trac.sagemath.org:sage into tmp</code>
</td></tr><tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=1294bf20bb6fe915a5371318cb8f9feacfa3f86a"><span class="icon"></span>1294bf2</a></td><td><code>trac #16964: review</code>
</td></tr></table>
TicketgitSat, 28 Feb 2015 14:50:25 GMTcommit changed
https://trac.sagemath.org/ticket/16964#comment:17
https://trac.sagemath.org/ticket/16964#comment:17
<ul>
<li><strong>commit</strong>
changed from <em>1294bf20bb6fe915a5371318cb8f9feacfa3f86a</em> to <em>e6a5a00acfb25935788ee67fc2defe28060343be</em>
</li>
</ul>
<p>
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
</p>
<table class="wiki">
<tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=e6a5a00acfb25935788ee67fc2defe28060343be"><span class="icon"></span>e6a5a00</a></td><td><code>trac #16964: review</code>
</td></tr></table>
TicketgitSat, 28 Feb 2015 15:15:29 GMTcommit changed
https://trac.sagemath.org/ticket/16964#comment:18
https://trac.sagemath.org/ticket/16964#comment:18
<ul>
<li><strong>commit</strong>
changed from <em>e6a5a00acfb25935788ee67fc2defe28060343be</em> to <em>0bf7dba36d3935c23ac385de7c28cf6d85055cd5</em>
</li>
</ul>
<p>
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
</p>
<table class="wiki">
<tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=8c0cf63c0e86d0cfc042149901ef9cfcde4d14ac"><span class="icon"></span>8c0cf63</a></td><td><code>trac #17863: remove stuff from src/ext</code>
</td></tr><tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=dfc820d0c69193066b0b42edb355e76ae84cbf91"><span class="icon"></span>dfc820d</a></td><td><code>trac #16964: merge sage-6.6.beta1</code>
</td></tr><tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=0bf7dba36d3935c23ac385de7c28cf6d85055cd5"><span class="icon"></span>0bf7dba</a></td><td><code>trac #16964: review</code>
</td></tr></table>
TicketgitSat, 28 Feb 2015 15:17:04 GMTcommit changed
https://trac.sagemath.org/ticket/16964#comment:19
https://trac.sagemath.org/ticket/16964#comment:19
<ul>
<li><strong>commit</strong>
changed from <em>0bf7dba36d3935c23ac385de7c28cf6d85055cd5</em> to <em>731ec907ee9cdc81bb06d647ebf0be26bd6ba117</em>
</li>
</ul>
<p>
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
</p>
<table class="wiki">
<tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=4f082b48627bc86c659487f4f84d1e36e37e5e3f"><span class="icon"></span>4f082b4</a></td><td><code>trac #16964: merge sage-6.6.beta1</code>
</td></tr><tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=731ec907ee9cdc81bb06d647ebf0be26bd6ba117"><span class="icon"></span>731ec90</a></td><td><code>trac #16964: review</code>
</td></tr></table>
TicketVincent DelecroixSat, 28 Feb 2015 15:17:35 GMT
https://trac.sagemath.org/ticket/16964#comment:20
https://trac.sagemath.org/ticket/16964#comment:20
<p>
Sorry I messed up the commits! Now everything looks good and doctest are fine...
</p>
TicketgitSat, 28 Feb 2015 17:25:50 GMTcommit changed
https://trac.sagemath.org/ticket/16964#comment:21
https://trac.sagemath.org/ticket/16964#comment:21
<ul>
<li><strong>commit</strong>
changed from <em>731ec907ee9cdc81bb06d647ebf0be26bd6ba117</em> to <em>297c68ad67657dc0e4c04a362827934a2e8ed17c</em>
</li>
</ul>
<p>
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
</p>
<table class="wiki">
<tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=a96f3f1def854014324ad94c4d108f9cd229b575"><span class="icon"></span>a96f3f1</a></td><td><code>trac #16964: review</code>
</td></tr><tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=297c68ad67657dc0e4c04a362827934a2e8ed17c"><span class="icon"></span>297c68a</a></td><td><code>trac #16964: (review) add an example</code>
</td></tr></table>
TicketMartin von GagernSat, 28 Feb 2015 18:40:21 GMTbranch changed
https://trac.sagemath.org/ticket/16964#comment:22
https://trac.sagemath.org/ticket/16964#comment:22
<ul>
<li><strong>branch</strong>
changed from <em>u/vdelecroix/16964</em> to <em>u/gagern/16964</em>
</li>
</ul>
TicketMartin von GagernSat, 28 Feb 2015 18:45:59 GMTcommit changed
https://trac.sagemath.org/ticket/16964#comment:23
https://trac.sagemath.org/ticket/16964#comment:23
<ul>
<li><strong>commit</strong>
changed from <em>297c68ad67657dc0e4c04a362827934a2e8ed17c</em> to <em>16f62a2ea8e36bf3ac9f781677e404e01366a594</em>
</li>
</ul>
<p>
I like how you avoid unneccessary minpoly computation. But I think we can do better than you did, by not having the intervals span negative and positive imaginary parts, but instead considering the absolute value of the imaginary part. I did something along these lines, but I now see that I'll have to rebase that on your latest forced push…
</p>
<hr />
<p>
New commits:
</p>
<table class="wiki">
<tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=731ec907ee9cdc81bb06d647ebf0be26bd6ba117"><span class="icon"></span>731ec90</a></td><td><code>trac #16964: review</code>
</td></tr><tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=16f62a2ea8e36bf3ac9f781677e404e01366a594"><span class="icon"></span>16f62a2</a></td><td><code>trac #16964: Use absolute value of imaginary part</code>
</td></tr></table>
TicketVincent DelecroixSat, 28 Feb 2015 18:52:56 GMT
https://trac.sagemath.org/ticket/16964#comment:24
https://trac.sagemath.org/ticket/16964#comment:24
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/16964#comment:23" title="Comment 23">gagern</a>:
</p>
<blockquote class="citation">
<p>
I like how you avoid unneccessary minpoly computation. But I think we can do better than you did, by not having the intervals span negative and positive imaginary parts, but instead considering the absolute value of the imaginary part. I did something along these lines, but I now see that I'll have to rebase that on your latest forced push…
</p>
</blockquote>
<p>
sorry for that... Your version is much simpler by the way!
</p>
<p>
But I am not completely convinced that this is optimal. Do you know if there is a cheap way to assert if two roots of a given (irreducible) polynomial have the same real value? In the example I added in the commit 297c68a <code>sorted(p2.roots(QQbar,False)</code> still takes hours.
</p>
<p>
Vincent
</p>
TicketVincent DelecroixSat, 28 Feb 2015 19:02:33 GMT
https://trac.sagemath.org/ticket/16964#comment:25
https://trac.sagemath.org/ticket/16964#comment:25
<p>
I can also rebase my changes on yours BTW.
</p>
TicketMartin von GagernSat, 28 Feb 2015 19:42:49 GMT
https://trac.sagemath.org/ticket/16964#comment:26
https://trac.sagemath.org/ticket/16964#comment:26
<p>
Would I be correct to assume that all this <code>ii_minus</code> and <code>ii_plus</code> handling in a96f3f1 is essentially equivalent to my use of absolute values in 16f62a2? If that is the case, then I'd rather drop a96f3f1 from the history, and instead rebase 297c68a onto 16f62a2. How does that sound to you?
</p>
TicketgitSat, 28 Feb 2015 19:46:17 GMTcommit changed
https://trac.sagemath.org/ticket/16964#comment:27
https://trac.sagemath.org/ticket/16964#comment:27
<ul>
<li><strong>commit</strong>
changed from <em>16f62a2ea8e36bf3ac9f781677e404e01366a594</em> to <em>aaf3eb648e1feca6e27c6dd3ffc5fdd9cee98191</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=aaf3eb648e1feca6e27c6dd3ffc5fdd9cee98191"><span class="icon"></span>aaf3eb6</a></td><td><code>trac #16964: (review) add an example</code>
</td></tr></table>
TicketVincent DelecroixSun, 01 Mar 2015 00:37:04 GMTcommit, branch changed; reviewer set
https://trac.sagemath.org/ticket/16964#comment:28
https://trac.sagemath.org/ticket/16964#comment:28
<ul>
<li><strong>commit</strong>
changed from <em>aaf3eb648e1feca6e27c6dd3ffc5fdd9cee98191</em> to <em>3f4afef46ab6a042cb2678394031cb5c26d89b1a</em>
</li>
<li><strong>branch</strong>
changed from <em>u/gagern/16964</em> to <em>u/vdelecroix/16964</em>
</li>
<li><strong>reviewer</strong>
set to <em>Vincent Delecroix</em>
</li>
</ul>
<p>
I choose a more meaningful test. Could you have a look?
</p>
<p>
As a consequence of your changes, we get great improvements on comparing rationals!!
</p>
<p>
new timings
</p>
<pre class="wiki">sage: a = QQbar(1/3)
sage: b = QQbar(1/2)
sage: %timeit cmp(a,b)
1000000 loops, best of 3: 881 ns per loop
</pre><p>
old timings
</p>
<pre class="wiki">sage: a = QQbar(1/3)
sage: b = QQbar(1/2)
sage: %timeit cmp(a,b)
10000 loops, best of 3: 34.4 µs per loop
</pre><p>
This is completely crazy.
</p>
<p>
If you are happy with my changes you can set to positive review.
</p>
<p>
Vincent
</p>
<hr />
<p>
New commits:
</p>
<table class="wiki">
<tr><td><a class="ext-link" href="http://git.sagemath.org/sage.git/commit/?id=3f4afef46ab6a042cb2678394031cb5c26d89b1a"><span class="icon"></span>3f4afef</a></td><td><code>trac #16964: better doctest</code>
</td></tr></table>
TicketMartin von GagernSun, 01 Mar 2015 00:58:04 GMTstatus changed
https://trac.sagemath.org/ticket/16964#comment:29
https://trac.sagemath.org/ticket/16964#comment:29
<ul>
<li><strong>status</strong>
changed from <em>needs_review</em> to <em>positive_review</em>
</li>
</ul>
<p>
I like it. Thanks for the review!
</p>
<p>
The speed gain for the rational numbers appear to be due to the fact that we no longer have to construct an algebraic number representation for the real part. <code>a._value.real()</code> is a lot faster than <code>a.real()</code>:
</p>
<pre class="wiki">sage: %timeit a._value.real()
1000000 loops, best of 3: 201 ns per loop
sage: %timeit a.real()
100000 loops, best of 3: 10.7 µs per loop
</pre><p>
(I still hope that someone, some day, will make all this work here obsolete by coming up with a better way to compare QQbar elements even if they don't share a minpoly. No idea how, though.)
</p>
TicketJeroen DemeyerSun, 01 Mar 2015 09:06:43 GMT
https://trac.sagemath.org/ticket/16964#comment:30
https://trac.sagemath.org/ticket/16964#comment:30
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/16964#comment:24" title="Comment 24">vdelecroix</a>:
</p>
<blockquote class="citation">
<p>
But I am not completely convinced that this is optimal. Do you know if there is a cheap way to assert if two roots of a given (irreducible) polynomial have the same real value?
</p>
</blockquote>
<p>
Using resultants, you can find a polynomial which has (a - b) as a root. And then you use interval arithmetic to find which root. Determining whether the imaginary parts are equal then is equivalent to checking whether a real polynomial has a pure imaginary root. The latter can probably be done using some Sturm-like algorithm.
</p>
TicketMartin von GagernSun, 01 Mar 2015 21:19:08 GMT
https://trac.sagemath.org/ticket/16964#comment:31
https://trac.sagemath.org/ticket/16964#comment:31
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/16964#comment:30" title="Comment 30">jdemeyer</a>:
</p>
<blockquote class="citation">
<p>
Using resultants, you can find a polynomial which has (a - b) as a root.
</p>
</blockquote>
<p>
Up to now I had assumed that most arithmetic in QQbar would eventually be performed using resultants. But it seems I was mistaken.
</p>
<pre class="wiki">sage: x = polygen(ZZ)
sage: p1 = x^5 + 6*x^4 - 42*x^3 - 142*x^2 + 467*x + 422
sage: p2 = p1((x-1)^2)
sage: r1 = QQbar.polynomial_root(p2, CIF(1, (2.1, 2.2)))
sage: r2 = QQbar.polynomial_root(p2, CIF(1, (2.8, 2.9)))
sage: a,b = polygens(QQ, 'a,b')
sage: %time p3 = r1.minpoly()(a + b).resultant(r2.minpoly()(b), b)
CPU times: user 62 ms, sys: 0 ns, total: 62 ms
Wall time: 68 ms
sage: [r for f in p3.factor()
....: for r in f[0].univariate_polynomial().roots(QQbar, False)
....: if r._value.overlaps(r1._value - r2._value)]
[-0.7266314748516305?*I]
sage: %time p4 = (r1 - r2).minpoly()
</pre><p>
One possible root of <code>p3</code> is <code>b=r2</code> and <code>a+b=r1</code> which means <code>a=r1-r2</code>. So eliminating <code>b</code> we get a (reducible, not minimal) polynomial in <code>a</code> which has that difference as one of its roots, just as you indicated. I try to identify that by looking at the roots <code>r</code> of the factors <code>f</code>, checking whether they overlap the numeric interval. The single result I obtain has zero real part, thus indicating that we should sort by imaginary part.
</p>
<p>
I lost patience waiting for that final result of <code>p4</code>, which should do pretty much the same thing in my opinion. My question is, why is it computed the way it is? Why do arithmetic operators for algebraic numbers compute some costly unions of number fields (which I believe is what they are doing), instead of using resultants to describe their results? And should we start some major rewrite effort to change that, i.e. to base most if not all arithmetic operations on resultants?
</p>
<p>
I can think of two possible problems. One is that we might be dealing with a special case here, and that perhaps number field unions are in general cheaper than resultants. If that were the case, what does that mean for speeding up comparisons? Another possible problem I can imagine is that the resultant could factor into several distinct polynomials, some of which might share a root. If that were the case, numeric refinement wouldn't be able to help choosing the right factor. Should we perhaps not factor the resultant polynomial, but instead compute roots for the fully expanded form?
</p>
<p>
I get the impression that this might trigger some major work. I hope that the reviewed changes can land as they are, while we investigate (in some other branch) how to tackle this more generic approach.
</p>
TicketMartin von GagernMon, 02 Mar 2015 17:06:13 GMT
https://trac.sagemath.org/ticket/16964#comment:32
https://trac.sagemath.org/ticket/16964#comment:32
<p>
Replying to <a class="ticket" href="https://trac.sagemath.org/ticket/16964#comment:31" title="Comment 31">gagern</a>:
</p>
<blockquote class="citation">
<p>
Why do arithmetic operators for algebraic numbers compute some costly unions of number fields (which I believe is what they are doing), instead of using resultants to describe their results? And should we start some major rewrite effort to change that, i.e. to base most if not all arithmetic operations on resultants?
</p>
<p>
I get the impression that this might trigger some major work. I hope that the reviewed changes can land as they are, while we investigate (in some other branch) how to tackle this more generic approach.
</p>
</blockquote>
<p>
Just created <a class="new ticket" href="https://trac.sagemath.org/ticket/17886" title="#17886: enhancement: Faster qqbar operations using resultants (new)">#17886</a> about using resultants to speed up most qqbar operations.
</p>
TicketVolker BraunTue, 03 Mar 2015 00:25:18 GMTstatus, branch changed; resolution set
https://trac.sagemath.org/ticket/16964#comment:33
https://trac.sagemath.org/ticket/16964#comment:33
<ul>
<li><strong>status</strong>
changed from <em>positive_review</em> to <em>closed</em>
</li>
<li><strong>resolution</strong>
set to <em>fixed</em>
</li>
<li><strong>branch</strong>
changed from <em>u/vdelecroix/16964</em> to <em>3f4afef46ab6a042cb2678394031cb5c26d89b1a</em>
</li>
</ul>
Ticket